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PREFACE 


This  final  report  represents  the  work  conducted  under  Contract  FA9453-13-C-0202  from  27 
December,  2012  to  22  July,  2016.  This  project  was  directed  by  Principal  Investigators  (Pis)  Dr. 
Kyle  T.  Alffiend  and  Dr.  Srinivas  R.  Vadali  of  Texas  A&M  University.  The  Pis  were  assisted  by 
Research  Associates  Dr.  Hui  Yan  and  Dr.  Kohei  Fujimoto,  and  PhD  students,  Bharat  Mahajan 
and  Lt  Col  Kirk  W.  Johnson. 

The  work  presented  in  this  report  is  a  summary  of  the  following  set  of  papers,  easily  accessible 
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1.  Yan,  H.,  Vadali,  S.R.,  and  Alffiend,  K.T.,  “State  Transition  Matrix  for  Relative  Motion 
Including  General  Gravitational  Perturbations,”  AAS  15-339,  AAS/AIAA  Space  Flight 
Mechanics  Meeting,  Williamsburg,  VA,  Jan  11-14,  2015. 

2.  Yan,  H.,  Vadali,  S.R.,  and  Alfriend,  K.T.,  “A  Recursive  Formulation  of  the  Satellite 
Perturbed  Relative  Motion  Problem,”  AAS/AIAA  Astrodynamics  Meeting,  San  Diego, 
CA,  Aug  2014. 

3.  Yan,  H.,  Vadali,  S.R.,  and  Alffiend,  K.T.,  “State  Transition  Matrix  for  Relative  Motion 
Including  J2  and  Third-Body  Perturbations,”  AAS  14-379,  AAS/AIAA  Space  Flight 
Mechanics  Meeting,  Santa  Fe,  NM,  January  2014. 

4.  Yan,  H.,  Vadali,  S.R.,  and  Alfriend,  K.T.,  “State  Transition  Matrix  for  Relative  Motion 
including  Higher-order  Gravity  Perturbations,”  AAS  13-793,  AAS/AIAA  Astrodynamics 
Meeting,  Hilton  Head,  SC,  Aug  12-15,  2013. 

5.  Mahajan,  B.,  Vadali,  S.  R.,  and  Alfriend,  K.  T.,  “Analytic  Solution  for  Satellite  Relative 
Motion  with  Zonal  Gravity  Perturbations,”  No.  15-705,  AAS/AIAA  Astrodynamics 
Specialist  Conference,  Vail,  CO,  August,  2015. 

6.  Mahajan,  B.,  Vadali,  S.  R.,  and  Alfriend,  K.  T.,  “Analytic  Solution  for  Satellite  Relative 
Motion:  The  Complete  Zonal  Gravitational  Problem,”  No.  16-262,  26th  AAS/AIAA 
Spaceflight  Mechanics  Meeting,  Napa,  California,  February,  2016. 

7.  Johnson,  K.  W.,  Vadali,  S.  R.,  and  Alffiend,  K.  T.,  "Comparison  of  Orbit  Element  Sets  for 
Modeling  Perturbed  Satellite  Relative  Motion,"  26th  AAS/AIAA  Space  Flight  Mechanics 
Meeting,  Napa,  California,  2016. 

8.  Fujimoto,  K.,  Alfriend,  K.  T.  and  Vadali,  S.  R.  “Bridging  Dynamical  Modeling  Effort  and 
Sensor  Accuracy  in  Relative  Spacecraft  Navigation,”  Presented  at  the  2015  AAS/AIAA 
Astrodynamics  Specialist  Conference,  Vail,  CO,  2015. 
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1  SUMMARY 


This  report  documents  the  work  performed  for  the  project  entitled  “Relative  Motion  Modeling  and 
Autonomous  Navigation  Accuracy”  being  performed  at  Texas  A&M  University  under  contract  to 
the  Air  Force  Research  Laboratory.  The  Pis  of  this  project  are  Drs.  Kyle  T.  Alfriend  and  Srinivas 
R.  Vadali. 

The  objectives  of  this  project  are: 

•  To  develop  a  methodology  for  determining  the  required  accuracy  of  the  dynamic  model 
based  on  system  requirements,  the  relative  navigation  accuracy  and  thruster  accuracy. 

•  To  extend  the  Gim-Alfriend  state  transition  matrix  (GA-STM)  for  relative  motion  to 
include  non-conservative  forces  and  the  effect  of  higher  order  gravitational  terms  on  the 
mean  motion. 

•  Validate  the  developed  models  using  numerical  simulation. 

In  current  practice,  the  dynamical  model  in  a  spacecraft  navigation  algorithm  is  often  set  ad  hoc 
without  explicit  regard  for  the  level  of  measurement,  guidance,  or  control  errors  expected.  The 
dynamic  model  for  the  relative  motion  of  two  satellites  should  be  consistent  with  the  accuracy  of 
the  relative  navigation  system  and  the  control  system,  otherwise  fuel  will  be  wasted  or  unnecessary 
computation  will  be  performed.  For  example,  including  more  perturbations  in  the  relative  motion 
dynamic  model  that  improve  the  predicted  motion  over  a  few  orbits  by  a  few  centimeters  provides 
no  real  improvement  if  the  relative  navigation  accuracy  is  no  better  than  10’s  of  centimeters.  The 
purpose  here  is  to  develop  a  methodology  that  would  simplify  the  workflow  of  designing 
navigation  systems  so  that  the  trade  space  between  navigation  system  parameters  and  dynamical 
model  fidelity  could  be  quickly  surveyed  in  lieu  of  performing  massive  numerical  simulations  for 
numerous  scenarios  and  system  parameter  variations.  First,  the  cost  of  a  particular  dynamical 
model  within  the  navigation  algorithm  was  expressed  as  a  combination  of  the  computation  cost 
and  the  maneuver  impulse,  Av,  cost  to  rectify  the  expected  state  error.  A  confidence  value 
accompanies  the  Av  cost  through  its  Mahalanobis  distance  with  respect  to  the  state  estimate 
uncertainty.  Then  two  methods  are  introduced  to  reduce  the  cost  computation  burden.  First,  the 
state  transition  matrices  of  different  dynamical  models  were  approximated  as  functions  of  the  time 
derivatives  of  the  kinematic  states.  Then,  a  simplified  form  of  the  linear  sensitivity  of  maximum  a 
posteriori  state  estimates  with  respect  to  errors  in  the  state  transition  matrix  was  derived.  The 
accuracy  of  these  approximations  is  deemed  sufficient  through  a  Monte  Carlo  simulation  for  a 
wide  range  of  formation  geometries  in  low  Earth  orbit. 

The  GA-STM  was  the  first  development  to  include  the  effect  of  the  absolute  and  differential  effects 
of  the  equatorial  bulge  term,  J2,  in  the  linearized  equations  for  the  relative  motion  of  satellites.  The 
second  objective  of  this  project  is  to  improve  the  dynamic  modeling  of  the  relative  motion  by 
adding  the  effects  of  higher  order  gravitational  perturbations,  lunar-solar  effects  perturbations,  and 
the  non-conservative  perturbations  to  the  GA-STM.  This  has  been  achieved.  The  accuracy  of  the 
expanded  GA-STM  was  evaluated  using  numerical  simulations  with  the  Goddard  General  Mission 
Analysis  Tool  (GMAT). 
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The  results  show  that 


•  Including  the  effects  of  the  higher  order  geopotential  short  period  effects  on  the  calculation 
of  the  relative  initial  conditions  reduces  the  secular  drift.  The  J\  effects  are  the  most 
dominant. 

•  Including  the  secular  and  long  period  effects  of  the  lunar  perturbations  on  the  relative 
motion  of  high  altitude  satellites,  e.g.,  geosynchronous,  is  important. 

•  The  use  of  an  exponential  density  model  for  the  relative,  not  necessarily  the  absolute 
motion  of  the  Chief  or  reference  satellite,  improves  the  accuracy  of  the  relative  motion. 

Other  developments  are 

•  The  evaluation  of  using  other  variables,  such  as  Hoots  variables,  for  the  relative  motion, 
and 

•  A  generalized  formulation  for  computing  the  effects  of  higher  order  geopotential  terms  on 
the  relative  motion. 

The  results  of  this  program  will  improve  the  precise  maintenance  of  satellite  formations  and 
provide  an  approach  for  selecting  the  appropriate  dynamic  model  used  for  the  type  of  relative 
motion  mission.  This  will  save  costs  in  the  analysis  and  formation  design,  and  in  the  development 
of  on-board  software  because  the  dynamic  model  will  be  consistent  with  the  required  performance 
and  relative  navigation  system. 
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2  INTRODUCTION 


Accurate  analytic  modeling  of  inter-satellite  relative  motion  is  indispensable  for  fast  and  reliable 
prediction,  fuel-efficient  formation  maintenance,  space  situational  awareness,  and  proximity 
operations.  In  contrast  to  numerical  propagation,  analytic  or  semi-analytic  solutions  not  only  effect 
faster  long-term  propagation  of  the  system  states,  but  also  provide  valuable  insight  for  mission 
analysis  and  parametric  studies.  One  of  the  main  challenges  underlying  the  development  of  a 
satellite  relative  motion  theory  for  low  and  medium  Earth  orbits  is  the  modeling  of  the  effect  of 
nonspherical  Earth  gravity  and  other  perturbations  such  as  differential  drag,  third-body  gravitation 
and  solar  radiation  pressure.  For  guidance,  control,  and  navigation,  the  analytical  model  in  the 
form  of  a  state  transition  matrix  (STM)  is  most  desirable.  In  addition,  the  sensing  or  navigation 
modules  of  spacecraft  use  estimation  theory  to  provide  estimates  of  states,  which  are  derived  from 
Kalman  filters.  The  complexity  of  the  dynamic  model  implemented  in  a  Kalman  filter  to  estimate 
the  relative  motion  of  two  satellites  should  be  consistent  with  the  accuracy  of  the  relative 
navigation  system  and  the  control  system,  otherwise  fuel  will  be  wasted  or  unnecessary 
computation  will  be  performed.  For  example,  including  more  perturbations  in  the  relative  motion 
dynamic  model  that  improve  the  predicted  motion  over  a  few  orbits  by  a  few  centimeters  provides 
no  real  improvement  if  the  relative  navigation  accuracy  is  no  better  than  10’s  of  centimeters. 

This  report  presents  work  performed  to  address  the  two  main  issues  raised  above. 

2.1  Objective 

The  objectives  of  this  project  are: 

1.  To  extend  GA-STM  [1]  for  relative  motion  to  include  non-conservative  forces  and  the 
effect  of  higher  order  gravitational  terms. 

2.  To  develop  a  methodology  for  determining  the  required  accuracy  of  the  dynamic  model 
based  on  system  requirements,  the  relative  navigation  accuracy  and  thruster  accuracy. 

3.  Validate  the  developed  models  using  the  GMAT  [2], 

2.2  Tasks 

The  tasks  undertaken  under  this  contract  are: 

1.  Dynamic  Model  for  Earth  gravitational  perturbation:  This  task  extends  the  GA-STM  to 
include  the  mean  and  periodic  effects  of  the  higher  order  Earth  gravitational  perturbations, 
i.e.,  the  zonal  and  tesseral  harmonics. 

2.  Dynamic  Model  Expansion  to  include  non-Earth  gravitational  perturbations:  Inclusion  of 
differential  drag,  luni-solar  third  body,  and  solar  radiation  pressure  perturbations. 

3.  A  method  to  determine  consistency  between  the  dynamic  model,  navigation  accuracy  and 
thruster  accuracy. 
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2.3  Organization  of  the  Report 


Section  3  of  the  report  presents  the  methods,  assumptions,  and  procedures  used.  Section  3. 1  of  the 
report  presents  the  research  undertaken  to  develop  an  extended  GA-STM  for  Earth  gravitational 
perturbations.  Section  3.2  of  the  report  includes  the  development  of  the  STM  for  non-Earth 
gravitational  effects.  The  analysis  and  the  development  of  a  method  to  determine  consistency 
between  the  dynamic  model,  navigation  accuracy  and  thruster  Accuracy  is  presented  in  Section 
3.3.  Section  4  presents  the  results  and  discussion.  The  report  concludes  with  recommendations  for 
future  research,  the  subject  of  Section  5. 
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3  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 

3.1  Dynamical  Model  for  Earth  Gravitational  Perturbations 

This  section  presents  the  development  of  relative  motion  theories  to  model  zonal  and  tesseral 
harmonics.  The  principal  tool  used  for  developing  analytical  satellite  theories  is  the  method  of 
averaging,  either  by  formal  perturbation  theories  [3,  4]  or  the  perturbation  approach  of  Kaula  [5]. 
Lie-series  based  approaches  can  be  used  to  determine  long  period  and  short  period  effects  that 
are  closed  form  in  eccentricity  (at  least  for  the  zonals)  and  the  algorithm  is  extendable  to 
compute  effects  up  to  any  order.  In  contrast,  Kaula’s  approach  produces  compact  expressions  for 
first-order  short  period  effects  due  to  the  zonals  and  tesserals  using  eccentricity  functions,  which 
involve  infinite-series  in  eccentricity.  These  methods  and  their  applications  are  discussed  briefly 
in  the  following  subsections. 

3.1.1  Principles  of  Averaging 

The  principles  of  averaging,  as  applied  to  satellite  theory,  are  illustrated  by  using  the  zonal 
harmonics  expansion  of  the  geopotential.  For  the  zonal  problem,  the  geopotential  in  terms  of 
spherical  harmonics  can  be  expressed  as 


r 


1  ^ n 

71=2 


Pn(sin(4>)) 


(1) 


The  symbols  r,  0,  Re,  p.,]n  and  Pn  represent  the  satellite  radius,  its  geocentric  latitude,  radius  of 
Earth,  gravitational  parameter,  zonal  geopotential  coefficients,  and  the  Legendre  polynomial  of 
the  nth  degree,  respectively.  The  corresponding  Hamiltonian  can  be  expressed  as  a  power  series 
in  J2  using  a  combination  of  Delaunay  and  classical  orbital  elements.  The  Delaunay  and  classical 

elements  are  related:  L  =  -J/jJa,  G  =  Z,v/l-e2 ,  H  =  G cos/,  /  =  M,  g  =  co,  and  h  =  Q  . 

3. 1.1.1  Deprit’s  Meth od 

The  Hamiltonian  for  the  complete  zonal  problem  can  be  written  as 

H  =  H0  +  J2HX  +  ^ H2  (2) 


where 
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Ho  =  -l£ 

H1  =  ^{t)3P2(sm(<j>)) 


H 2  =  £~  3 


h2. 


=  En=3  7*Jnffi  (ar)U+1  Pn(sm(<t>)) 


The  above  Hamiltonian  H  for  the  complete  zonal  problem  can  be  normalized  using  Deprit’s  Lie- 
transform  based  canonical  perturbation  method.  This  approach  works  by  separating  the  secular 
and  periodic  effects  with  the  help  of  two  successive  Delaunay  normalizations  of  the  gravitational 
potential.  The  first  normalization  averages  out  the  short  period  terms  containing  mean  anomaly, 
l,  and  produces  a  short  period  generating  function.  The  resulting  singly-averaged  Hamiltonian 
includes  long  period  as  well  as  secular  terms.  The  second  normalization  is  necessary  to  separate 
the  long  period  terms  containing  the  argument  of  perigee  angle,  g,  to  produce  a  doubly-averaged 
Hamiltonian  (or  Kamiltonian)  consisting  of  only  the  secular  terms.  The  Lie-transform  based 
perturbation  equations  up  to  third-order  used  for  separating  the  secular  and  periodic  effects  can 
be  written  as 


H0  =  K0 
(H0,  Wi)  +  Hi  =  Ki 
(H0,  W2)  +  (JTi  +  Ku  W{)  +  H2  =  K2 
(Ho,  Wo)  +  (2ifi  +  Ku  W2)  +  (H2  +  2K2,  Wx)  -  ((K^Wx),  Wx)  +  H3  =  K3 


where  K0,  K1,  K2,  and  K3  are  the  Kamiltonians  of  successive  orders  and  similarly,  W1,W2,  and 
W3  are  the  generating  functions  of  successive  orders.  The  parentheses  represent  the  Poisson 
brackets,  which  are  evaluated  in  terms  of  the  canonical  Delaunay  element  set  D  = 

[l,  g,  h,  L,  G,  H].  The  total  Kamiltonian,  K,  up  to  order  3  and  the  generating  function,  W,  up  to 
order  2  are  respectively, 


K  =  K0  +  J2K\  +  31k2  +  ^K3  +  0(4) 

2!  3!  (4) 

72 

W  =  Wx  +  J2W2  +  -^W3  +  0(4) 

Because  J2  is  considered  as  the  only  first-order  perturbation,  the  first-order  Kamiltonian  K1  and 
short-period  generating  function  obtained  during  the  first  normalization  of  H,  have  no 
contributions  from  the  other  zonals.  The  second-order  perturbations  in  H2  appear  for  the  first 
time  at  the  second  order  in  the  perturbation  equations  and  involves  computation  of  Poisson 
brackets  only  at  the  third-order.  As  a  result,  the  second-order  contributions  K2  n  to  the  total 
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Kamiltonian  K2  and  W2  n  to  the  total  short-period  generating  function  W2  due  to  any  zonal 
harmonic  Jn  (n  >  3)  can  be  computed  as 


W2,n  = 


(5) 


The  subscript  p  of  the  integral  sign  indicates  that  the  integration  is  performed  for  the  periodic 
part  of  the  integrand  only,  i.e.,  terms  involving  cosine  and  sine  functions  of  true  or  mean 
anomaly.  By  making  use  of  definition  of  the  perturbing  potential,  the  above  integrals  can  be 
written  as 


1  2! 

^2,n  — 


/J.R 


27T  J|  an+1\/l  —  e2  Jo 


/  n  \  n—1 

t?  L  (;)  p"(sin^))  d/ 


w2,n  =  t/^  ( K2,M  +  ^ 


Jfan+W T 


n  f  /n\n~ 1 

J  (;)  p»(sin(«)  v 


(«) 


The  above  integrals  can  be  computed  in  closed  form  for  an  arbitrary  zonal  to  produce 
expressions  for  the  single-averaged  Hamiltonian  and  short  period  generating  function.  A  similar 
procedure  can  be  used  to  perform  the  second  normalization  to  average  out  the  long-period  terms 
dependent  on  the  perigee  angle,  g,  to  derive  the  long  period  generating  function.  Once  the 
expressions  for  averaged  Hamiltonian  and  generating  functions  are  available,  the  secular  as  well 
as  periodic  effects  on  a  satellite  due  to  these  perturbations  can  be  easily  computed.  The  secular 
variations  of  the  orbital  elements  can  be  derived  using  Hamilton’s  canonical  equations,  which 
propagates  the  mean  elements,  £m  in  time  using  the  singly  or  doubly-averaged  Hamiltonian.  The 
periodic  contributions  to  the  mean  elements  up  to  second-order  can  be  computed  using  the 
following  near-identity  transformation  to  get  the  osculating  elements,  £ 0 : 

t2 

So  =  Sm  +/2(£m,  Wi)  +ji(((£m'W1)tW1)  +  (£m,W2j)  +  0(7!)  (?) 

3. 1.1.2  Kaula’s  method 

Kaula’s  approach  [4]  treats  the  full  geopotential  expansion,  including  zonal  and  tesseral  harmonics. 
The  potential  is  expressed  in  terms  of  the  classical  elements  as 

u  =  ;  +  sr=2  EU.  Sp=o  w  0  elM  O )slmpq  (<u.  m,  ft,  eCST )  (8) 
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The  S  functions  contain  the  gravitational  coefficients  and  the  F  and  G  functions  are  the  inclination 
and  eccentricity  functions,  respectively,  and 

^  _  Qm  cos  0  +  Sim  sin  0  <  if  (/  —  m)  even 

—Sim  cos  0  +  CLm  sin  0 ,  if  (/  —  m)  odd  ' 


where 


0  =  ®impq  =  (!-  2p)o>  +  (l  -  2p  +  q)M  +  m(H  -  6GST ) 

where  /  is  the  zonal  degree,  m  is  the  tesseral  order,  p  is  the  index  of  inclination  function  F,  q 
is  the  power  of  the  eccentricity  function  G  and  0GST  is  the  Greenwich  Sidereal  Time. 

Explicit  forms  of  the  F  and  G  functions  have  been  tabulated  in  Kaula  [5].  The  averaged 
Hamiltonian  for  Ji  is  obtained  by  setting  in  Eq.  (8),  l  =  2,  m  =  0,  p  =  1,  q  =  0  and  the 
averaged  Hamiltonian  for  J4  is  the  potential  when  /  =  4,  m  =  0,  p  =  2,  q  =  0. 

Kaula’s  linear  solutions  are  useful  for  modeling  short  period  corrections  to  the  classical  orbital 
elements,  especially  for  orbits  of  small  eccentricity.  The  series  expansions  in  the  G  functions  can 
be  truncated  at  a  suitable  order  for  small  eccentricities.  According  Kaula’s  theory  [5],  the  first 
approximations  to  the  perturbations  of  the  orbital  elements  or  the  solutions  of  the  Lagrange 
planetary  equations  for  the  geopotential  model  are  expressed  as 


^ elmpq  ~  nai+3eQplmp^lpq[vO  2p  +  q)  (l 

hilmpq  =  nai+ 3??  sin  (g,  Flmp^lpq  [0  ~  2p)  COS  i  — 

A/,,  —  ^  r]  ly-1  r  d_^}vq  t  •  dFimt 

A0)lmpq  nal+3QVle  Flmp  gg  COt  1  di 

AO  —  d£lmp  r  <r 

m LlmPq  nal+3q  sin  i©  di  “iP^lmpq 

kMlmpq  =  naj+3^  [~ve  1  +  2(1  +  1  )Gipq] 


(10) 

(11) 

(12) 

(13) 

(14) 

(15) 


For  near  circular  orbits,  Kaula’s  solution  to  the  J2  problem  is  obtained  by  setting  q=0  and  l-2p+q± 
0: 
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Aasp  =  Aa2ooo  +  ^a2020  : 

*  3/2R|sln2'cos21 

2  a 

(16) 

Aisp  =  Ai2 ooo  +  ^2020  ~ 

3J2Re  sin  2i 

-  , —  cos  2A 

8a2 

(17) 

Afi5p  =  Afi2000  +  ^^2020 

=  3hRlfs,sm2A 

4  a2 

(18) 

3.1.2  The  Gim-Alfriend  State  Transition  Matrix 

An  analytic  satellite  theory  can  be  used  to  propagate  multiple  satellites  using  the  results  from  the 
previous  section.  Relative  motion  between  satellites  can  be  discerned  by  differencing  the  orbital 
elements  or  the  states  of  any  two  satellites  in  a  formation.  Even  though  this  approach  is 
reasonable  and  applicable  to  arbitrary  formation  sizes,  the  drawback  is  the  computational  time 
and  the  information  exchange  required.  In  addition,  a  direct  analytic  solution  for  the  relative 
motion  problem  also  benefits  in  formulating  guidance  and  control  algorithms  for  rendezvous  and 
proximity  operations.  For  formation  control,  especially  for  proximity  operations,  the  differential 
orbital  elements  of  a  deputy  can  be  represented  as  a  first-order  variation  of  the  orbital  elements 
of  the  reference  or  chief  satellite.  These  differential  orbital  elements  can  be  directly  propagated 
in  time  in  the  presence  of  the  gravitational  perturbations  using  the  differential  perturbation 
effects.  This  approach  originally  used  for  computing  the  state  transition  matrix  for  the  relative 
motion  in  the  presence  of  J2  perturbation  up  to  first-order,  is  known  as  GA-STM.  It  directly 
propagates  in  time  the  relative  states  of  the  deputy  in  a  curvilinear  frame  attached  to  the  chief. 

The  GA-STM  is  composed  of  three  different  matrices:  the  geometric  transformation  matrix  X, 
the  relative  mean- to-osculating  transformation  matrix  D,  and  the  relative  mean  STM  cp.  The 
relative  mean  STM  propagates  in  time  the  relative  mean  elements  of  the  deputy  with  respect  to 
the  chief,  using  the  mean  rates  of  the  chief’s  orbital  elements.  The  relative  mean  orbital  elements 
are  then  converted  into  the  relative  osculating  elements  when  operated  on  by  D.  Finally,  the 
geometric  transformation  matrix  transforms  the  relative  osculating  orbital  elements  into  the 
curvilinear  orbit  frame.  Using  these  three  matrices,  GA-STM,  ®,  can  be  written  as 

®(t,  t0)  =  X(t)£>(t)0(t,  to)D-1(to)X-1(t0)  (19) 


where  the  propagation  of  the  relative  position  and  velocity  states  in  the  curvilinear  coordinate 
system  is  achieved  from  time  t0  to  t. 

3. 1.2.1  The  Geometric  Transformation  Matrix 

The  geometric  transformation  matrix  converts  the  relative  orbital  elements  into  the  relative 
position  and  velocity  states  in  the  curvilinear  orbit  frame.  The  curvilinear  frame  is  defined  with 
respect  to  an  imaginary  circle  with  its  origin  at  the  chief  and  its  radius  taken  as  the  geocentric 
radius  of  the  chief.  The  relative  position  with  respect  to  the  chief  can  be  defined  in  terms  of  three 
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coordinates:  x  represents  the  difference  in  the  radii  of  the  two  satellites,  and  y  and  z  represent  the 
curvilinear  distances  from  the  chief,  along  and  normal  to  the  imaginary  circle,  respectively.  If  x 
represents  the  relative  states  in  the  curvilinear  frame  and  Ae  represents  the  relative  orbital 
elements,  then  the  geometric  transformation  matrix  is  defined  as 

x(t)  =  X(t)Ae(t)  (20) 

where  x  =  [x,  y,  z,  x,  y,  z]T  and  Ae  =  ed  —  ec.  The  subscripts  d  and  c  indicate  the  deputy  and 
chief,  respectively,  and  the  boldface  letters  represent  vectors.  The  geometric  transformation 
matrix  can  be  computed  by  taking  the  first  variation  of  the  chief’s  radial  distance  and  its 
direction  cosine  matrix.  The  first  three  rows  of  E  represent  the  position  equations  and  are  given 
as  follows  [1]: 


M  c  =  [SRc]c  +  Cci8Cji[Rc]c  (21) 

where  [Rc]c  =  [ftc  0  0]r,Mc  =  [x  y  z]T ,  and  Cci  =  R3(f  +  a))R1(i')R3(Sl')  . 

The  notation  [  ]c  represents  a  vector  with  components  expressed  in  the  chief’s  orbital  frame  and 
the  prefix  6  represents  the  variation  of  the  quantity.  The  symbol  R  represents  the  simple  rotation 
matrix,  with  the  subscript  specifying  the  type  of  rotation.  In  a  similar  fashion,  the  relative 
velocity  of  the  deputy  in  the  curvilinear  frame  can  also  be  computed  by  making  use  of  the 
transport  equation: 


[x  +  mxx]c  =  [dta  +  CciSCTci[Vc]c  (22) 

where  [Vc]c  =  [Vr  Vt  Vn]T .  The  symbols,  Vr,  Vt,  and  Vn  represent  the  chief’s  velocity 
components  in  the  radial,  tangential  and  normal  directions,  respectively.  The  symbol  m 
represents  the  angular  velocity  of  the  chief  satellite,  which  can  be  expressed  in  terms  of  the 
osculating  orbital  elements  as 


u 7  = 


flsini 
sin(/  +  co) 


(23) 


1  d  R 

where  fl  = - and  Rn  =  I? Hi. 

nab  sin  i  di  P  J  L  1 

The  symbols,  p,  n,  and  b  represent  the  parameter  of  the  chief’s  orbit,  mean  motion,  and  the 
semiminor  axis,  respectively.  The  symbol  Rp  represents  the  perturbing  potential  and  is  needed  to 
compute  the  nodal  angle  rate.  In  the  case  of  original  GA-STM,  Rp  consists  of  J2  effects  only, 
however  the  effects  due  to  an  arbitrary  zonal  harmonic  can  be  incorporated  by  including  the 
corresponding  perturbing  potential  in  Rp.  The  equations  for  relative  position  x  and  relative 
velocity  x  can  be  expressed  in  terms  of  the  orbital  elements.  Expressions  for  the  elements  of  X  in 
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terms  of  nonsingular  and  equinoctial  elements  for  the  first-order  J2  problem  are  given  in  [1],  The 
inverse  of  I  is  required  to  transform  the  initial  relative  osculating  states  in  the  curvilinear  frame 
to  the  relative  osculating  orbital  elements;  it  can  either  be  computed  numerically  or  by  reversion, 
using  symbolic  algebra. 

3. 1.2. 2  The  Relative  Mean-to-Osculating  Transformation  Matrix 

The  matrix  D  converts  the  relative  mean  elements  into  the  relative  osculating  elements.  It 
captures  the  short-period  and  long-period  effects  on  the  differential  orbital  elements  and  is 
defined  as  the  Jacobian  of  the  contact  transformations  as  given  in  Eq.  (7).  The  D  matrix  can  be 
computed  using  the  following  relations 


D  —  DspDlp 


(24) 


where  =  [^]  and  D„  =  gg]. 

In  the  above  relations,  the  symbols  £m,  £0,  £lP  represent  the  chief’s  mean  elements,  the 
osculating  elements  and  the  elements  with  short-period  effects  averaged  out,  respectively. 

3.1. 2.3  The  Relative  Mean  STM 

The  relative  mean  STM ,(p,  propagates  in  time  the  relative  mean  differential  elements  of  the 
deputy  relative  to  those  of  the  chief.  It  models  the  secular  effects  on  the  differential  mean 
elements  by  using  the  first-order  variation  of  the  mean  rates  of  the  chief.  The  following  equations 
define  the  relative  mean  STM: 


d£m(t)  =  4>(t,t0)d£m(t0)  (25) 

where  <p(t,  t0)  =  a£ni(t)  and  the  symbols  £m(t)  and  £m(t0)  represent  the  chief’s  mean  elements 

d£m(to) 

at  time  t  and  initial  time  t0,  respectively. 

3.1.3  Extended  Relative  Motion  Model  based  on  Nonsingular  Elements 

The  addition  of  the  higher-order  corrections  to  the  orbit  theory  allows  one  to  extend  the 
capabilities  of  the  GA-STM.  Brouwer’s  theory  [3]  incorporates  the  secular,  short  periodic 
long  periodic  effects ./?,  and  the  secular  and  long  periodic  effects  of  Jr  and  J3-J6.  For  the 
geopotential,  the  magnitudes  of  J22  and  the  higher  coefficients  are  of  the  same  order.  Thes 
effects  have  been  used  to  extend  the  GA-STM  significantly.  In  addition,  the  short  period  1 
of  J2 2  have  also  been  incorporated.  The  nonsingular  elements  are:  [a,  A  =  M  +  a),  i,  q1  = 
e  cos(o)),  q2  =  e  sin(co)  ,h  =  D], 
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3. 1.3.1  Short  Period  Terms  due  to  J22 

Kinoshita  [6],  in  his  third-order  theory,  provides  an  expression  for  the  short-periodic  generating 
function  of  0(  J21 )  valid  for  small  eccentricity.  His  expressions  for  the  mean-osculating 
transformation  can  be  simplified  considerably  by  setting  eccentricity  to  zero  and  retaining 
second-order  accuracy.  These  correction  terms  for  the  nonsingular  elements  are 


«f22=F[(“i3s,'2+^si4)cos5;l+ 

9  173  _  449  ,  81  75  ,  231  , 

(—  —  +  — —  si2  —  — —  si4)  cos  3A  +  (—  — —si2  +  — —  si4)  cos  A] 
v  16  64  128  J  v16  8  64  J  J 


sp22  J2  r 

h  =7fi[ 


sp22  J2  r(  21  51  .4)  . 

^  =I5[l-64SI  +128s'  jSln5A  + 

9  155  ,  395  ,  63  21  ,  465  „ 

(“  16  +  ~MSi  ~  128 Si  )  Sin  31  +  (16  “  TSi  +'MSi  )  Sinl] 


,,  /!  (9  3  39  A  /  15  ,  99  x 

1  =  ZJ 1 [{32  ~  64*  ~  64Si  )  $in 4A  +  (~TSi  +  32s1*)  Sin  2 


,,  /|  /  9  3  ,\  9  , 

hsp22  =  73-  ci  [(  —  —  —  —  si2  )  sin  4A  +  —  si2  sin  2A1 

L8  LV  32  32  /  8  J 

LSP22  =  Zl  [—si4  cos4A  +  (—si2  —  —  si4)  cos  2A 
L7  64  V  8  8  / 

(9  45  ,  141  A 

+  (s-i6s'2  +  lTsi4)l 


where  si  =  sin  i,  ci  =  cos  i.  From  these  expressions,  the  corrections  to  the  other  elements  are 

asp22  =  2  Lsp22Va 


;SP 22  _ 


^sp  22 


Asp22  +  2 qs^22  sin  A  —  cos  A 
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3. 1.3.2  STM  with  Kaula ’s  Theory  for  Zonal  and  Tesseral  Harmonics 

Assuming  small  eccentricity,  the  Kaula’s  corrections  for  the  nonsingular  elements  are  computed 
as  [7]: 


A almpq  2  nai+2Q  FlmpGlpq (J  R^lmpq 

[mmpd-%* (jfj)  + 

2(1  +  l)F[mpGlpq] 

^ hmpq  nal+3q  sin  i©  ^lmp^lpq  [(^  2p)  COS  i  (32) 


^Rllmpq  ^^Impq  COSO)  (<?A(0  ^impq  Sin  Ct) 

^Cfalmpq  ^C[mpq  sin  CO  +  (<?A 0)fmpq  COS  CO 

AO  tlRe  dflmv  F 

Impq  nal+3q  sin  i©  di  U^lmpq 

For  the  numerical  example,  gravity  coefficients  of  the  20x20  gravity  field  are  obtained  from  the 
JGM-2  Model.  The  above  equations  are  used  to  obtain  the  orbital  elements  of  the  chief  and 
deputy,  using  their  respective  mean  elements.  The  unit  sphere  approach  is  used  to  compute  the 
relative  motion  variables  in  the  local  vertical  local  horizontal  frame.  Equivalent  results  from 
GMAT  are  also  obtained  for  the  20x20  gravity  field. 

3.1.4  STM  for  the  Complete  Zonal  Perturbation  Problem 

The  nonspherical  gravitational  problem,  including  all  the  zonal  harmonics  is  referred  to  as  the 
complete  zonal  problem.  This  subsection  describes  in  brief  the  methodology  for  deriving  the 
STM  for  satellite  relative  motion  that  includes  the  perturbation  effects,  closed-form  in 
eccentricity,  due  to  zonal  harmonics  up  to  an  arbitrary  degree.  The  secular  as  well  as  the  periodic 
effects  up  to  second-order  of  any  zonal  harmonic  on  the  orbital  elements  of  a  satellite  are 
modeled  using  the  Deprit’s  method.  By  making  use  of  the  GA-STM  framework  described  in  the 
previous  subsection,  these  effects  were  then  incorporated  into  a  STM  solution  for  the  relative 
motion  problem.  Secular  or  mean  rates  up  to  3rd  order,  long-period  and  short-period  effects  up  to 
2nd  order  for  two  different  sets  of  orbital  elements  have  been  computed:  the  nonsingular  set 
(nonsingular  for  zero  eccentricity)  and  the  equinoctial  set.  The  nonsingular  set  is  singular  for 
equatorial  orbits,  but  the  equinoctial  element  set  is  completely  nonsingular  for  equatorial  and 
circular  orbits.  It  is  noted  that  the  expressions  for  secular  and  periodic  variations  of  the 
equinoctial  elements  are  significantly  larger  than  those  for  the  nonsingular  set. 

A  complete  second-order  analytic  solution  for  an  artificial  satellite  as  well  as  a  STM  for  relative 
motion  incorporating  the  perturbation  effects  due  to  J2  —  /6  has  been  computed  for  the 
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nonsingular  element  set  [9].  Even  for  nonsingular  elements,  the  expressions  for  secular  and 
periodic  effects,  especially  short-period  effects  at  2nd  order,  are  too  large  and  cumbersome  for 
computations  by  hand.  Therefore,  the  Maple  symbolic  algebra  package  was  used  to  compute  the 
STM  in  terms  of  nonsingular  element  sets.  Because  of  the  use  of  nonsingular  elements,  this 
version  of  relative  motion  STM  has  singularities  for  reference  orbits  that  lie  in  the  equatorial 
plane.  Additionally,  extending  this  STM  beyond  J6  by  computing  expressions  explicitly  for  each 
zonal  proved  to  be  a  very  time  consuming  and  computer  resource  expensive  process,  even  with 
the  use  of  a  symbolic  algebra  package. 

To  address  these  limitations,  a  slightly  different  approach  was  used  to  compute  the  STM  in  terms 
of  equinoctial  elements.  This  approach  extended  the  earlier  work  done  by  Saedeleer  [8],  in  which 
generalized  expressions  for  first-order  averaged  Hamiltonian  and  generating  function  for  the 
short-period  effects,  closed  form  in  eccentricity,  due  to  an  arbitrary  zonal  harmonic  were  given. 
These  formulae  were  originally  derived  for  the  classical  orbital  elements.  Using  the  Deprit’s 
method,  generalized  analytic  formulae  were  derived  for  second-order  secular  and  short-period 
and  first-order  long-period  variations  of  the  equinoctial  elements  due  to  an  arbitrary  zonal 
harmonic  Jn  (n  >  3).  To  derive  the  STM,  analytic  formulae  for  partial  derivatives  of  the  mean 
rates,  short-period  and  long-period  transformations  with  respect  to  equinoctial  elements  were 
also  derived  by  hand.  All  these  formulae  are  valid  for  any  elliptic  reference  orbit  without  any 
singularities  related  to  zero  eccentricity  or  inclination  values.  Because  J2  causes  the  dominant 
effect,  it  is  considered  as  a  first-order  perturbation  while  all  the  higher  zonals  as  second-order 
perturbations.  To  validate  the  accuracy  of  the  proposed  STM,  results  were  compared  with  direct 
numerical  propagation  using  GMAT.  The  following  subsection  provides  the  analytic  expressions 
for  propagating  a  single  satellite  used  to  propagate  the  equinoctial  elements  of  the  chief.  It  is 
noted  that  the  analytic  propagation  of  the  chief  is  needed  to  compute  the  STM  for  the  relative 
motion. 


3. 1. 4. 1  Secular  Effects 

This  subsection  presents  the  secular  effects  on  the  equinoctial  elements  of  an  artificial  satellite 
for  the  complete  zonal  problem.  For  any  harmonic  Jn  (n  >  3),  the  generalized  analytic  formulae 
for  second-order  mean  rates  are  provided  here.  For  J2,  the  expressions  for  secular  rates  up  to 
third-order  are  explicitly  computed  using  Maple.  By  using  the  analytic  formulae,  the 
contributions  to  the  secular  effects  due  to  the  zonal  harmonics  Jn  (n  >  3)  can  be  conveniently 
added  to  the  J2  secular  rates  without  having  to  go  through  the  process  of  two  Delaunay 
normalizations  of  the  perturbing  potential.  In  order  to  derive  the  secular  effects  due  to  an 
arbitrary  zonal,  the  closed- form  expression  for  the  doubly-averaged  Hamiltonian  or  Kamiltonian 
K2  n  is  sought.  Equation  (6)  provides  the  integrals  necessary  to  be  evaluated  for  finding 
expressions  for  K2n.  Appendix  A.  1  provides  helpful  formulae  for  evaluating  the  integrals  in  Eq. 
(6). 

Using  the  expressions  for  K2  n  ,  the  formulae  for  the  secular  or  mean  rates  of  the  orbital  elements 
can  be  derived  by  using  the  Hamilton’s  canonical  equations.  It  is  noted  that  the  Equinoctial 
elements  are  not  canonical.  Therefore,  the  secular  rates  of  the  Delaunay  elements  are  computed 
first,  which  can  then  be  used  to  compute  the  rates  of  the  Equinoctial  elements.  Using  the 
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canonical  equations  of  motion,  the  formulae  for  the  mean  rates  of  Delaunay  elements  D  due  to 
even  harmonics  are  computed  as  follows: 
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LfJ 
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The  total  secular  rates  of  the  Delaunay  elements,  D,  are  found  by  adding  the  secular  rates  due  to 
J2,  T)j2,  to  the  contributions  of  the  higher  degree  zonal  harmonics,  Dn,from  the  above  formulae, 
as  shown  below: 


t2  00 

i>  =  vj2  +  -£  V  t>„  (34) 

n— 3 

The  Equinoctial  elements  (a ,  e,  and  i  represent  the  three  classical  orbital  elements:  semimajor 
axis,  eccentricity  and  inclination,  respectively.)  are 
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a  =  a 


A  —  l  g  h 
Pi 


tan  (  -  J  cos (h) 


p2  =  tan  J  sin (h) 

qi  =  e  cos  (9  +  h ) 
q2  =  e  sin  (9  +  h ) 


(35) 


The  Equinoctial  elements  of  the  reference  satellite  can  be  propagated  from  the  initial  time  t0  to 
time  t  as  follows: 


a(t)  =  a(t0 ) 

A  (t)  =  A(to)  +  (i ^  +  9  +  h)(t  —  to) 

Pi{t)  =  Pi  (to)  cos  (h(t  -  t0))  -  p2{to)  sin  (h(t  -  t0)) 

P2(t)  =  P2(to)  cos (h(t  -  to))  +  Pi  (to)  sin(A(t  -  t0)) 

Qi(t)  =  9i (to)  cos((p  +  h)(t-  t0))  -  92 (to)  sin((p  +  h)(t  - 10)) 

92(t)  =  92 (to)  cos((p  +  h)(t-  t0))  +  9i (to)  sin ((9  +  h)(t  -  t0)) 


3. 1.4.2  Short  Period  Effects 

Using  the  expansion  formulae  given  in  Appendix  A.  1 ,  the  analytic  expression  for  the  second- 
order  short-period  generating  function  W2iTl,  closed  form  in  eccentricity,  corresponding  to  an 
arbitrary  zonal  Jn  (n  >  3)  can  be  computed  from  Eq.  (6).  This  expression  for  W2  n  can  be  written 
succinctly  as: 


W2, 


n 


=  Wi  (a,  e,  i,  /,  /,  g)  +  W2(a,  e,  i)  W3(/,  9) 


(37) 
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In  the  above  equations,  the  square-bracketed  terms  applies  for  odd  harmonics  and  c  =  n  —  2j  — 
2  p.  The  formulae  for  second-order  short  period  contributions  due  to  any  zonal  harmonic  ]n  (n  > 
3)  can  now  be  derived  by  evaluating  the  Poisson  brackets  with  equinoctial  elements.  The 
formulae  for  these  short-period  contributions  are  derived  by  hand  in  this  work  and  are  given  as 
below: 


(a,  W2,n)  —  —2 


+  W2 


dj_ 

dl 


(43) 
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^  i  .  a  iff  1  \zsLGH 
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In  the  above  equations,  frequently  occurring  expressions  are  replaced  by  defining  new  terms  to 
save  space.  The  definition  of  these  terms  are  provided  in  Appendix  A.2.  The  expressions  for 
short-period  corrections  up  to  second-order  due  to  J2  are  computed  explicitly  by  using  the  Maple 
software.  The  osculating  equinoctial  elements  £0  can  be  computed  by  using  the  following  near¬ 
identity  transformations: 


eo  =  eLP  +  J2  (SLp,  Wf  p)  +  £  (((£LP,  ttfp)  ,  Hf p)  +  (Slp,  Wfp ))  +  0(4) 

(49) 

Slp  =  S0-  -h  {So,  wfp)  +  §  (((So,  wfp ) ,  Wfp)  -  (So,  W2SP))  +  0(4) 

In  the  above  equations,  Wfp  represents  the  first-order  generating  function  due  to  J2  and  W2P 
represents  the  complete  second-order  generating  function  due  to  J2  and  higher  zonals.  The 
second  of  the  above  equations  is  the  inverse  near-identity  transformation. 

3. 1.4.3  Long  Period  Effects 

Once  short-period  effects  are  averaged  out,  a  second  Delaunay  normalization  of  the  singly- 
averaged  Hamiltonian  is  required  to  compute  the  long-period  generating  function.  Since  the 
long-period  generating  function  is  not  a  function  of  l,  the  perturbation  equations  given  in  Eq.  (3) 
simplify  because  the  first  Poisson  bracket  involving  H0  vanishes.  As  a  result,  the  first-order  long- 
period  generating  function  W£p  can  be  computed  using  the  following  second-order  perturbation 
equation 


(H1+K1,W[p)+H2  =  K2 


(50) 


In  the  above  equation,  H2  represents  the  second-order  singly-averaged  Hamiltonian  that  includes 
long-period  terms  dependent  on  g  due  to  the  zonals  ]n  (n  >  3),  K2  is  the  second-order 
Kamiltonian  with  only  the  secular  terms.  Because  the  J2  potential  has  no  long-period  terms 
dependent  on  g,H1=K1.  The  Poisson  bracket  in  the  above  equation  can  be  evaluated  to  result  in 
the  following  equation  for  Wppx  for  an  arbitrary  zonal  hamronic  Jn  (n  >  3). 

f  K2,n  d9  (51) 

z  dG  JlP 

The  subscript  Ip  on  the  integral  sign  denotes  the  integration  of  the  long-periodic  terms  of  the 
integrand  only.  The  first-order  Kamiltonian  Kx  has  contributions  from  J2  only  and  its  value  is 
given  as  [9]. 


1  Re2p  (3  cos2  (i)  —  l) 
4  a3?;3 


(52) 
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Substituting  the  above  value  for  K1  and  the  formulae  for  K2n,  the  formulae  for  the  first-order 
long-period  generating  function  for  even  zonal  harmonics  (n  >  3)  can  be  written  as 


wi,n  =  sin(i)Wi(a,e,i)W2(5) 


(53) 


where 
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(57) 


Using  the  above  expression  for  the  long-period  generating  function,  the  analytic  formulae  for  the 
long-period  effects  due  to  any  zonal  harmonic  Jn  (n  >  3)  can  be  computed  by  evaluating  the 
Poisson  brackets  for  the  equinoctial  elements.  These  formulae  for  the  first-order  long-period 
effects  are  derived  and  given  as  follows: 


(a,  =  0 


(58) 
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LP\  .  ...  Wi(a,e,z) 
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It  is  noted  that  the  order  of  terms  in  a  product  in  the  above  equations  must  be  preserved  because 
of  presence  of  the  summation  indices.  The  long-period  effects  can  be  added  to  the  mean 
Equinoctial  elements  £m  to  compute  Equinoctial  elements  with  short-periods  effects  averaged 
out,  SLP,  using  the  following  near- identity  transformation  up  to  second-order. 


Slp  =  £m  +  J2  (£m,  Wfp)  +  (((£m,  w£p)  ,  w£p)  +  (£m,  WPP))  +  0(j|) 

2  (64) 

£,«  =  £lp  -  h  (£lp,  Wtp)  +  %  (((Slp,  W,lp) ,  Wtp)  -  (£lp,  W2lp))  +  0(4) 
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3. 1.4.4  STM  for  the  Perturbed  Relative  Motion 

The  perturbations  effects  computed  in  the  previous  subsections  due  to  an  arbitrary  zonal 
harmonic  Jn  (n  >  3)  were  incorporated  into  the  GA-STM  framework  by  computing  the  partial 
derivatives  of  the  formulae  for  mean  rates,  long-period  and  short-period  effects  with  respect  to 
the  equinoctial  elements.  These  analytic  formulae  for  the  partial  derivatives  are  derived  by  hand 
and  implemented  in  a  MATLAB  code,  however  they  not  included  in  this  report  to  save  space. 
The  partial  derivatives  of  the  mean  rates  are  incorporated  into  the  differential  mean  STM, 
whereas  the  partial  derivatives  of  the  long-period  and  short-period  formulae  are  used  to  update 
the  Dlp  and  DSP  matrices  of  the  GA-STM  as  defined  in  the  previous  subsections.  To  update  the 
geometric  transformation  matrix  given  in  Subsection  3. 1.2.1,  the  chief’s  nodal  angle  rate 
formulae  must  be  updated  to  include  the  effects  of  higher  zonals  by  using  the  following 
perturbing  potential  (This  equation  is  not  numbered.) 


A  =  Mi +2^2  (65) 

3.1.5  Gim-Alfriend  State  Transition  Matrix  in  Terms  of  Hoots  Elements 

This  portion  of  the  study  reformulated  the  GA-STM  for  linearized  satellite  relative  motion  in 
terms  of  Hoots  elements  [10].  The  Hoots  elements  are  y1  =  r,y2=  f,  y3  =  rf,  y4  = 
sin  i/2  sin  u,  y5  =  sin  i/2  cos  u,  and  y6  =  u  +  h,  where  r  is  the  orbit  radius,  /  is  the  true 
anomaly,  i  is  the  inclination,  h  is  the  right  ascension  of  the  ascending  node,  u  =  f  +  g  is  the 
argument  of  latitude,  and  g  is  the  argument  of  perigee.  The  Hoots  elements  are  completely 
nonsingular  for  all  eccentricities  and  inclinations. 

3. 1. 5. 1  Sensitivity  Matrix  D 

The  sensitivity  matrix  D  (t)  is  defined  as  D  =  dy/dy",  the  Jacobian  of  the  osculating  elements 

with  respect  to  the  mean  elements.  Considering  the  osculating  elements  as  a  vector  of  nonlinear 
functions  of  the  mean  elements,  and  expanding  the  elements  for  a  deputy  satellite  as  a  Taylor 
series  about  the  reference  or  chief  satellite,  yields  the  following: 


Td(0  =  yc(t)  + 


dy 


dy' 


(y"a(0-y"c(0) 


+ 


y"c 


(66) 


This  shows  that  D  can  function  as  a  linearized  operator  for  relative  motion,  mapping  the 
differences  in  mean  Hoots  elements  into  differential  osculating  Hoots  elements: 


£y(t)  =  D(t)[y"a(t)-y"c(0] 


(67) 
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Importantly,  this  relationship  can  be  inverted  at  the  initial  time,  providing  a  way  to  convert  the 
osculating  initial  conditions  to  mean  relative  initial  conditions: 


sy’(to)  =  D  1(t0)  \yd(t  o)  -  yc(t  o)]  (68) 

Given  correction  terms  y  —  y"  from  the  mean-to-osculating  transformation  of  a  particular 

d(y-yrr) 

perturbation  theory,  we  can  compute  the  sensitivity  matrix  as  D  =  /6x6  +  v~  ~  ,  where  I  is  the 

identity  matrix.  If  the  perturbation  theory  uses  successive  transformations  (for  example,  long- 
period  and  short-period  corrections),  then  we  can  apply  the  D  operators  in  succession: 


D  =  DspDlp  = 


dy  dyi 
dyr  dy" 


*6X6 


d(y-y')N 


dyr 


*6X6 


+ 


d(y'-y")N 


dyr 


(69) 


It  is  possible  to  compute  D  based  on  the  first-order  /2-only  corrections  from  Hoots  theory 
(equivalent  to  Brouwer  theory,  but  with  singularities  removed  and  certain  other  improvements) 

as  D  =  /6x6  +  +  D^sp\  where  Z)0p)  =  -,  D(sp:>  =  and  the  product 

j)(sp)[)(ip)  |ias  been  neglected  as  second-order  in  J2.  Once  the  correction  terms  y'  —  y"  and  y  — 
y'  are  expressed  in  terms  of  only  the  Hoots  elements,  forming  matrices  Z)6p)  and  D^sp^  is  a 
straightforward  matter  of  taking  partial  derivatives.  Note  that  using  mean  elements  y",  rather 
than  intermediate  elements  y',  as  inputs  when  computing  D(sp^  would  constitute  yet  another 
approximation,  introducing  error  of  second  order  in  J2. 

3. 1.5. 2  Mean  State  Transition  Matrix  <p~ 

Considering  the  state  y"(t )  as  a  vector  function  of  the  initial  conditions,  and  expanding  the  state 
for  a  deputy  satellite  as  a  Taylor  series 

yVO  =  y"c(t) 

or  8y"(t) 


about  the  chief,  yields  the  following: 


dy'Xt ) 


df(t0) 

dy"(t) 


y"c(t  0) 


(y"i(to)-y"c(to))  + 


dy"(t0) 


8y"(t0)  + 


(70) 

(71) 


y"c(fo ) 
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is  the  linearized  relative  state  transition  matrix 


This  shows  that  <p-(t,  t0) 


dy"(t) 


dy"(t0 ) 


y"c(t  o) 


(STM),  so  that  (neglecting  terms  of  second  order  in  the  relative  coordinates)  8y"{t )  = 

<p-(t,  t0)5y"(t0).  Finding  the  partial  derivatives  of  y"(t)  with  respect  to  the  initial  conditions  is 
most  conveniently  accomplished  through  a  change  of  variables  to  a  set  of  elements  x:  x1  =  a, 
x2  =  e  sin  l,  x3  =  e  cos  l,  x4  =  sin  i/2  sin  h,  x5  =  sin  i/2  cos  h,  and  x6  =  l  +  g  +  h,  where  e  is 
eccentricity  and  l  is  mean  anomaly.  Then  we  can  say  that 


dy"(t)  _  dy"(t)  dx"(t)  dx"(t0) 
dfiuj  ~  dx"(t)  dx"(t0)  dy"(t0) 


Hoots  found  the  partial  derivatives  of  y  with  respect  to  x  ,  and  these  can  be  assembled  into  J  = 


the  geometric  transformation  between  the  two  sets,  which  has  the  same  form  whether 

dyrr  dyr  dy 

computed  using  mean,  intermediate,  or  osculating  elements  (i.e.,  =  -^).  The  vector  x 

is  expressed  in  terms  of  the  initial  conditions  as 


x/(t)  =  *2  fro)  cos[(n"  +  lp)(t  -  t0)]  +  X3  (t0)  sin [{n"  +  Z/)(t  -  t0)] 

*3(0  =  -*2  (fo)  sin[(n"  +  i/)(t  -  t0)]  +  x'/ (t0)  cos[(n"  +  i/)(t  -  t0)] 

(73) 

< (t)  =  *4  (t0)  cos[/r/(t  -  t0)]  +  x'5'  (t0)  sin[hp'(t  -  t0)] 
x's(t)  =  -x4  (t0)  sin[A/(t  -  t0)]  +  X5  (t0)  cos[Ap'(t  -  t0)] 
x6”(t)  =  x6"(t0)  +  (n"  +  ip  +  gp  +  hp)(t  -  t0) 
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where  n  =  /p/a3  is  the  mean  motion  (so  that  n"  =  /p/x/'^o)3)  and  the  subscript  p  indicates 
the  secular  rates  due  to  perturbations  (which  can  be  found  from  Brouwer  theory  [3]).  Let 


<Px(t’  to) 


dy_"(t) 
dy"(t0 ) 


.  Then 


<Py(t,to)  =  J(t)<p-(t,  t0)J  1(t0)  (74) 

where  all  initial  conditions  are  for  the  chief  satellite.  Note  that  x  and  <p-(t,  t0)  are  nonlinear 

functions  of  the  perturbed  mean  secular  rates  and  will  have  to  be  re-derived  to  account  for  J22  or 
higher  zonal  harmonics. 


3. 1.5.3  Relative  Transformation  Map  I 


Reference  [1]  reports  expressions  for  the  spherical  curvilinear  coordinates  of  a  deputy  satellite’s 
relative  position  r  (t)  and  velocity  v(t)  in  terms  of  the  chief’s  osculating  nonsingular  elements 
(a,  u,  i,  q1  =  e  cos  g,  q2  =  e  sin  g,  and  h )  and  the  deputy’s  relative  osculating  nonsingular 
elements  (8a,  8u,  Si,  8q±,  Sq2,  and  8h).  These  can  be  easily  mapped  into  expressions  in  terms 
of  chief  and  relative  classical  orbit  elements  using  the  variations,  8q1  =  cos  g  8e  —  e  sin  g  8g, 

8q2  =  sin g  8e  +  e  cos  g  8g,  8u  =  8g  +  8f,  and  8f  =  ^8e  +  ^ 81 .  Some  terms  in  v(t) 

depend  on  the  perturbed  rates  of  change  in  u  and  h,  which  can  be  found  from  Gauss’s  Variational 
Equations  in  terms  of  the  perturbing  acceleration  vector  due  to  the  Earth’s  gravitational  zonal 
harmonics.  The  expressions  for  r  (t)  and  v(t)  can  then  be  transformed  into  Hoots  elements  via 
the  following  relationships  (all  derived  from  the  definitions  of  the  Hoots  elements  given  above): 


sin  it  = 


V4 


COS  u  = 


V5 


,  sin /  =  ™  s/  =  i (2!L-  l\  a  =  sini/2  = 

nae  e  \  y±  J 


sini/25  ~  sini/25 

Vt42  +  T52»  and  cos  i/2  =  /  l  —  y42  —  y52,  where  g  =  Vl  —  e2.  Solving  these  relations  for  the 
final  variable,  e,  yields  a  4th-order  polynomial  whose  only  root  on  e  6  [0,1)  is  e  = 

T.Jy  i2y34  +  yi2y22y32  —  2py1y32  +  p2.  The  differential  elements  are  transformed  using  a 

linearized  map  containing  the  partial  derivatives  of  the  classical  orbit  elements  with  respect  to 
the  Hoots  elements.  Finally,  the  coefficients  of  the  differential  Hoots  elements  are  formed  into 
matrix  T(t),  so  that 


fc(t)r  £(t)T  =  KOW)  (75) 


The  portion  of  the  transformation  due  to  the  perturbing  acceleration  can  be  partitioned  into  a 
separate  map  B(t),  so  that  T(t)  =  A(t)  +  B(t ). 
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3.2  Dynamical  Model  Expansion:  Non-Earth  Gravitational  Perturbations 

The  dominant  perturbation  on  satellite  motion  for  objects  in  LEO  is  the  equatorial  bulge  term,  J2. 
Atmospheric  drag  affects  objects  up  to  about  1000  km  altitude.  At  higher  altitudes,  particularly, 
geosynchronous  altitude,  the  perturbations  resulting  due  to  the  moon  and  sun  can  be  as  large  as 
the  Earth  gravitational  effects  and  need  to  be  evaluated.  Also,  at  these  altitudes  the  solar  radiation 
effects,  which  are  not  a  function  of  altitude,  can  also  have  a  significant  effect.  At  altitudes  under 
1000  km  atmospheric  drag  effects  can  be  significant  and  need  to  be  considered.  The  effect  of 
these  perturbations  on  the  relative  emotion  of  satellites  are  addressed  in  this  chapter.  Lunar  and 
solar  effects  are  addressed  in  Section  3.2.1,  solar  radiation  effects  are  covered  in  Section  3.2.2 
and  atmospheric  drag  in  Section  3.2.3.  These  effects  are  incorporated  into  the  Gim-Alfriend 
STM. 

3.2.1  Dynamic  Model  Expansion:  Third  Body  Perturbations 

The  effects  of  lunar  perturbations  on  satellites  have  been  studied  extensively,  using  perturbation 
methods  and  averaging.  The  first  lunisolar  disturbing  function  was  developed  for  the  secular  and 
long  periodic  terms  in  the  1950s  by  Kozai  [11]  and  expanded  by  Musen  et  al  [12].  Kaula  [5] 
investigated  a  general  method  to  represent  the  disturbing  potentials  in  terms  of  orbital  elements. 
The  advantage  of  this  approach  is  that  specific  terms  can  be  studied  in  a  convenient  and  general 
manner.  For  example,  for  the  secular  motion,  the  averaged  Hamiltonian  is  the  potential  obtained 
by  setting  the  coefficients  depending  on  the  satellite  mean  anomaly  and  argument  of  perigee  to 
zero.  The  method  was  revisited  by  Giacaglia  [13]  and  a  treatment  is  presented  in  Vallado  [14]. 
Although  this  infinite  series  summation  method  provides  a  general  formulation,  efficient 
numerical  implementations  require  recursive  formulations  [15].  Kozai  [11]  developed  an 
alternative  method  to  compute  the  lunisolar  perturbations.  The  disturbing  function  was  expressed 
in  terms  of  the  orbital  elements  of  the  satellite  and  the  polar  geocentric  coordinates  of  the  Sun 
and  the  Moon.  From  the  disturbing  function  form,  the  short  periodic  terms  can  be  eliminated  by 
taking  the  average  with  respect  to  the  mean  anomaly  of  satellite  motion.  The  summation  of  the 
terms  in  the  potential  is  carried  up  to  order  five.  Prado  [16]  applied  a  double  averaging 
technique  to  obtain  the  third-body  disturbing  potential  using  the  formulation  of  the  restricted 
three-body  problem.  Since  the  short  periodic  terms  are  removed  from  the  model,  it  is  convenient 
for  the  study  of  long-term  orbit  stability  and  leads  to  fast  computations.  Recently,  the  accuracy  of 
the  doubly-averaged  model  has  been  improved  by  including  the  lunar  orbit’s  eccentricity  and 
inclination  in  the  studies  presented  in  References  [17-19]. 

Since  the  perturbed  relative  motion  problem  is  very  complex,  our  effort  has  focused  on 
simplifying  the  dynamic  models  of  relative  motion.  This  section  considers  the  third-body 
perturbation  to  extend  the  fidelity  of  the  GA  STM,  based  on  the  previous  work  by  Roscoe,  Vadali 
and  Alfriend  [12,  21],  The  second-order  effect  of  the  averaged  lunar  perturbation  is  modeled  with 
nonsingular  elements.  For  nonlinear  simulations,  the  relative  motion  of  the  Moon  with  respect  to 
the  Earth  is  incorporated  from  ephemerides  data  obtained  from  the  Jet  Propulsion  Laboratory’s 
HORIZONS  website.  The  initial  conditions  for  the  averaged  model  are  determined  by  a  least 
squares  method  developed  in  Reference  [18].  The  nonlinear  solution  is  validated  by  the  General 
Mission  Analysis  Tool  (GMAT)  [2], 
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3. 2. 1. 1  Averaged  Disturbing  Potential 

In  the  following,  the  subscript  “m”  refers  to  the  variable  related  to  the  Moon  and  the  symbols  “S” 
and  “C”,  written  with  subscripts,  stand  for  sine  and  cosine  functions,  respectively.  The  disturbing 
potential  due  to  the  Moon  is 


PrnG(me+mm) 

sjr2  +rl-2rrm  cos  5' 


(76) 


where  pim  =  mm/ ( me  +  mm),  G  is  the  gravitational  constant,  me  is  the  mass  of  Earth  mass,  >nm 
is  the  mass  of  the  Moon,  r  is  the  distance  from  Earth  to  a  satellite,  and  rm  is  the  Earth-Moon 
distance.  The  included  angle  between  the  vectors  r  and  rm  ,  denoted  by  S ,  can  be  obtained  from 
the  dot  product  of  the  unit  vectors  r  and  rm. 

In  the  Earth  centered  inertial  (ECI)  frame,  the  unit  vector  rm  can  be  expressed  in  terms  of  the 
components 


Gem  ^nm  Cim 

ym=SnCdm+CnmSeCim 

K  =SA 


The  ECI  components  can  be  expressed  in  the  perifocal  frame  as 


x  =(CnC  -SnS  C.)x  +(SnC  +CnS  C.)y  +S  S.z 

mp  \  Q  co  0.  0)1  /  m  \  Q  co  Q  co  i  )  J  m  co  i  m 

y  =(-CaS  -SaS  C.)x  +(-SaS  +CnC  C.)y  +C  S.z 

y  mp  V  O  oj  O  co  i  )  m  \  O  oj  Q  qj  i  /  y  m  oj  i  m 

z  =  SaSx  -CaSy  +C-Z 

TYl  12  l  TYl  12  l  171  l  171 


In  the  perifocal  frame,  the  unit  vector  r  has  components 


xn  =  cos  f 

P  J 

yP  =sin  / 

4=° 


where  /  is  the  true  anomaly.  From  Eqs.  (77-79),  we  have 

cos  S  =  a  cos  /  +  P  sin  / 


(77) 


(78) 


(79) 


(80) 


where 


«  =  (Cm Ca>  ~  SMiSoA)Cem  +  (CAa^CCm  +  SM1C0JChn  +  S(0SiSim)Sdn 
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where  M  is  the  mean  anomaly  and  F  is  an  arbitrary  potential  function.  Taking  the  first  average 
over  the  mean  anomaly  of  the  satellite  orbit  [7]  yields  the  averaged  potential 


The  identities  (So  )  =  (C%  )  =  -  and  (So  Ca  )  =  0  have  been  used  to  derive  the  above 
equation. 

The  second  average  is  taken  with  respect  to  the  lunar  period.  Considering  the  lunar  orbit’s 

3 

eccentricity,  we  have  from  [20],  ((— )3)  =  (1  —  e^)-2.  Hence,  the  second  average  of  Eq.  (83)  is 

rm 

((R2  ))  =  Mm"”a  (l  -  e2m  p  [(6  +  9e2 )  A  + 1 5 e2B  -  12e2  -  8]  (84) 

where 
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A  =  C2n(  1  +  CfCfJ  +  S2a(C£2  +  CfJ  +  SfSfm  +  i52£m52£CAa 


B  = 


Cla{  1  -  C£2Q2m)  +  ^n(Q2m  -  Q2)  -  SfSfm  - -S2imS2iCAn 


c 


2  oo 


—  C$2An Ci^im  ~  ^Aa^i^2im)^2co 


3.2.1. 2  Secular  and  Long  Period  Rates 

Lagrange’s  planetary  equations  can  be  obtained  for  the  double-averaged  potential  as 
a  =  0 

V  d((R2)) 


e  = 


na2e  do 


nazr]  sin  i  L  do  d£l  J 

1  d((R2)) 


n  =  — 1 

na2r]smi  di 

^  _  V  d((R2))  _  cos  i  d((R2)) 
na2e  de  na2r]s\ni  di 

M  =  n  t  9«R2)>  2  d««2» 

na2e  de  na  5a 

where  n  is  the  mean  motion  and  rj  =  Vl  —  e2.  The  required  partial  derivatives  are 

d«52» 


da  8 

d((52))  _  3/tmn^a2e 

de  8 


[(6  +  9e2)4  +  15e25  -  12e2  -  8] 


[321  +  55-4] 


di  16 

d«52»  _  BwPmOL2  r 

da  ”  16 

d«52»  _  ^mn^a2  r 

dco 
d2l 


16 


,  7X  d4  _  d5i 

(6  +  9e2)  —  +  15e2  — 
di  di  J 

,  dA  "  ,  dB - 

«  +  9^lm  +  15e  sa 

(6  +  9e2)|^+15e2^- 

0(1)  0(1). 


—  (^im  ^An  Qm  ^An)^2i  +  ^2im^2iQn 

d4  1 

^  =  (C4  “  Ci2  C4  “  ^i2)^2An  “  2  S2imS2iSAn 


(85) 


(86) 

(87) 

(88) 

(89) 

(90) 

(91) 

(92) 
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dB 
<9  00 


~  T  ^Atl^2i  S2iSim  •?2imC2iC^^]C2(0 

+(S2A  nsisim  +  SAnciS2im)S2(Ji 

$2A n(Q  ^im~^im  ~  Q  )  T  2  ^2im^2i^ASl  ^2u> 

— (2C2AaCiSim  —  C^a^i^2im)^2u> 

CL(  1  -  C?tfJ  +  S|„(C4  -  C,2)  -  S}Sl  -  is2i„,S2fC4n 

—2(S2AnCiSim  ~  S^a^i^2im)^2u> 


dA 

doo  _  ° 


an 


-2 


->20) 


(93) 

(94) 

(95) 

(96) 


The  doubly-averaged  lunar  perturbation  model  is  obtained  by  substituting  Eq.  (84)  into  the 
Lagrange’s  planetary  equations  (85). 

3.2.1. 3  Extended  GA  STM  Including  Third  Body  Perturbations 

The  development  of  the  STM  for  long  period  effects  of  the  3rd  body  perturbations  is  now 
addressed,  hence,  set  D  (t)  as  the  identity  matrix. 

3.2. 1 .3. 1  Secular  and  Long  Period  Terms 

The  STM  of  the  current  set  of  orbit  elements  with  respect  to  the  initial  orbit  elements, 

0e,  considering  only  the  secular  and  long  period  terms,  can  be  calculated  from  the  following 
relationships: 


A  —  A0  +  AAt 
i  =  /(j  +  /At 


( 


q[  =  e cos  to  =  1 

V 


15  H  n2  dB  A  A 
- —2! — —T\ At 

16  n  d(Q  j 


[qw  cos  A O)  -  q20  sin  Am) 


f 

q2  =  e  sin  ft)  =  1 


15  p  n2  dB  A 

16  n  d(0 


(g10  sin  Aft)  +  q20  cos  Aft)) 


Q  —  Qq  +  QA/ 


(97) 


where  the  subscript  “0”  means  the  initial  state  and  Aco  =  d)(t  —  t0).  The  STM  is  obtained  by 
assuming  that  the  orbit  element  rates  f l,  (6  ,  and  A  are  constant  and  that  A o>  is  small.  The  details 
of  the  derivations  are  presented  in  Appendix  B. 
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3.2.1 .3.2  Geometric  Transformation  Matrix 


The  relative  motion  is  best  visualized  in  a  local  vertical  local  horizontal  (LVLH)  frame  centered 
on  the  chief  or  a  curvilinear  system  [2],  The  I(t)  matrix  is  used  to  transform  the  osculating 
elements  to  LVLH  coordinates.  The  normal  acceleration  is 


h  sin  i  ^ 

1  i 

rsin# 


From  Lagrange  equations 


Q  = 


1  cffc, 

na2r\  sin  i  di 


(98) 


(99) 


3. 2.1. 4  Nonlinear  Simulation  Model  and  Initial  Conditions 

The  ECI  reference  epoch  is  set  at  the  J2000  epoch.  The  XY  plane  is  the  plane  of  the  Earth’s  orbit 
at  the  reference  epoch.  The  X  axis  points  to  the  ascending  node  of  the  instantaneous  plane  of  the 
Earth’s  orbit  and  Earth’s  mean  equator  at  the  reference  epoch.  The  Z  axis  is  perpendicular  to  the 
XY  plane  in  the  directional  sense  of  the  Earth’s  North  Pole  at  the  reference  epoch.  The  Y  axis  is 
determined  by  the  right-hand  rule.  The  simulations  include  the  gravitational  perturbations  and  the 
lunar  perturbations.  The  high  fidelity  ephemerides  data  for  the  lunar  motion  are  obtained  from 
the  Jet  Propulsion  Laboratory’s  HORIZONS  website;  the  data  is  then  interpolated  using  the 
method  described  in  Reference  23. 

Lagrange’s  planetary  equations  provide  the  rates  of  change  of  the  classical  orbital  elements  due 
to  the  doubly-averaged  lunar  perturbation.  However,  because  the  disturbing  potential  has  been 
averaged,  these  now  represent  the  rates  of  a  new  set  of  “lunar”-averaged  orbital  elements,  rather 
than  the  instantaneous  osculating  elements.  If  equation  (85)  is  applied  without  correcting  the 
initial  conditions  for  these  differences,  the  results  will  become  increasingly  inaccurate  as  the 
equations  are  propagated  forward  in  time.  In  the  absence  of  an  accurate  conversion  between  the 
osculating  and  averaged  elements,  a  least  squares  method  is  used  to  correct  the  initial  conditions. 

In  our  problem 


Ax  —  (aosc0  Q-is>@oscO  @lsXoscO  hs’QloscO  RllsO’  tf2osc0  R2ls0>^osc0 
—  Q-lso) 

y  =  (,ais>  9is>  hs>  <h  Is’  <i2is>  nto) 

f  Cv)  (floSC’^OSC’^OSC’  Qlosc>  92osc»^osc) 


(100) 

(101) 

(102) 


where  the  subscript  “osc”  stands  for  osculating  elements  and  “/.v”  represents  long  and  secular 
averaged  elements,  while  “0”  refers  to  initial  values. 
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3.2.2  Dynamic  Model  Expansion:  Solar  Radiation  Effects 

Solar  radiation  pressure  (SRP)  is  a  nonconservative  perturbation.  For  high  altitude  orbits,  solar 
radiation  pressure  can  be  the  dominant  perturbative  force,  particularly  for  satellites  with  large 
area  to  mass  ratio.  By  intentionally  aligning  a  satellite’s  orientation  with  respect  to  the  sun,  this 
perturbative  force  can  be  used  to  maneuver  the  vehicle.  For  many  years,  the  focus  of  SRP  studies 
was  confined  to  the  minimization  of  its  effects  on  spacecraft.  More  recently,  considerable  effort 
has  been  put  into  the  development  of  applications  to  exploit  the  SRP  effects  for  purposes  of 
interplanetary  propulsion,  namely  with  solar  sails  [22],  Zeng,  et  al.  [23]  presented  a  time-optimal 
trajectory  design  for  a  novel  dual-satellite  sailcraft  to  accomplish  mid  or  far-term  interstellar 
missions 

Application  of  SRP  for  formation  flying  has  been  investigated  by  various  researchers.  Williams 
and  Wang  [24]  considered  a  satellite  formation  with  a  solar  wing  and  it  was  shown  that  a  solar 
torque  can  be  generated  roughly  along  the  orbit  direction.  This  torque  can  prevent  the  secular 
out-of-plane  growth  in  a  low-Earth  orbit  formation  that  is  caused  by  differential  nodal  drift. 
Kumar,  et  al.  [25]  and  Gong,  et  al.  [26]  demonstrated  the  feasibility  of  using  SRP  for  maintaining 
the  desired  satellite  formation  for  different  scenarios. 

In  this  section  the  GA  STM  is  extended  to  include  the  differential  SRP  for  relative  motion.  The 
SRP  model  is  developed  and  the  short  and  long  period  dynamics  due  to  SRP  are  presented.  Only 
the  SRP  perturbation  is  considered.  The  differential  SRP  is  primarily  caused  by  the  differential 
area-to  mass  ratio  (AMR)  perturbations.  The  contribution  here  is  that  a  new  variable  AMR  is 
introduced  to  the  GA  STM  to  deal  with  the  AMR  perturbations.  Numerical  results  show  relative 
position  errors  are  significantly  reduced  after  the  differential  SRP  is  included  in  the  GA  STM. 

3. 2. 2. 1  Solar  Radiation  Pressure  Model 


The  acceleration  due  to  SRP  is  [14,  27] 


™cosH2[^COs4+(1-^ 


(103) 


where  P  is  the  solar  radiation  pressure  (P=4.56e-6  N/m2),  A  is  the  cross  sectional  area,  m  is 
mass,  (j)inc  is  the  incident  angle,  n  is  the  unit  vector  of  the  cross  sectional  plate,  s  is  the  unit 

vector  from  the  satellite  to  the  Sun,  and  Ca ,  Cd,  Cs  are  the  coefficients  of  absorption,  diffuse  and 
specular  reflectivity,  respectively.  Notice  Ca  +  Cd  +  Cs  =  1  and  assume  the  satellite  is  spherical 


a  = - ks  =  PkCj 


(104) 


where  k  =  1  +  Cs  + — Cd  and  the  area-mass-ratio  is  CB 


A_ 

m 
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The  components  of  the  acceleration  in  the  radial,  transverse  and  normal  direction  are 


Ur  =  -PkCB  [ur  cos  e  +  U,  sin  0) 

Ut  =  -PkCB  ( Ut  cos  6  -  Ur  sin  0}  (105) 

Uh=-PkCBUh 


U r  =  cos  Q  cos  Ae  +  sin  Q  sin  Ae  cos  s 

Ut=-  sin  Q  cos  i  cos  Ae  +  cos  Q  cos  i  sin  Ae  cos  s  +  sin  i  sin  Ae  sin  s 
Uh  =  sin  Q  sin  i  cos  Ae  -  cos  Q  sin  i  sin  Ae  cos  s  +  cos  i  sin  Ae  sin  s 


3.2.2.2  Secular,  Long  Period  Rates  and  Short  Period  Terms 
3.2.2.2.1  Secular  and  Long  Period  Rates 

Substituting  Eq.  (106)  into  the  nonsingular  Gauss  equation  yields 


da 

dt 

d6_ 

dt 

di 

dt 


2a~ 


h 


(<?!  sin  0-q2  cos 9)Ur+— Ut 

r 


h  r  sin  0  cos/ 
r2  h  sin  i 
r  cos# 


Uu 


h 


-Uu 


dq j 

dt 

dq2 

dt 


—\  sinOU,,  + 
h 

P_ 
h 


cos  0  H — ql 


LI 

PJ 

p  J 

- 

7  \ 

T 

cos  6Ur  + 

1+- 

sin<9H — q2 

l  p) 

P 

Ut+—q2  sin  0  cot  i  Uh 
P 


Ut - qx  sin  6 cot  iUh 

P 


dO. 

dt 

dCB 

dt 


_  r  sin  6 
h  sin  i  1 

=  0 


(106) 


(107) 
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Assuming  small  eccentricity  gives 

6  =  /l  +  2qx  sin  A  -  2q2  cos  A 
r 

—  =  1  -qx  cos  A  -q2  sin  A 
a 

sin  0  =  — q2  +  sin  A  +  qx  sin  2  A  -  q2  cos  2  A 
cos  0  =  —qx  +  cos  A  +  qx  cos  2  A  +  q2  sin  2 A 

Substituting  Eq.  (108)  into  Eq.  (107)  gives 


(108) 


a  =  -  PkCB  [Ur  sin  A  -  E7  cos  A  +  (Urqt  -  UqJ  sin  2  A  -  (Urq2  +  £/<?,)  cos  2A} 
i  =  — PkCBUh^3ql  -2cosA  -  qx  cos2A  -q2  sin2Aj 


2na 

Q  = - - - PkCU ,  {3 q  -2sin X  —  q.  sin2 A  + cos2A| 

2«asin/  5  *  L  J 

[3-1-  1-  .1/- 

q  = - -s  —U  +  — £/  cos2A  +  — £/  sm2A  +  —  ([/  q  +  3(7 g  jsinA 

[2  '  2  ?  2  r  4Vr  f 

\  _  _  3  _  _  3__  1 

+  —  [5Urq2  -3Lf  ^JcosAh-— ~~ ^rgJsin3A  +  — +  ^<Z1)cos3A  [  +  Qg2  cos/ 


3  -  U 


U 


q  = - —  <!  — U  -4 — Lsin2An — -  cos  2  A  +  —(3/7  g  -  5U  q  )  sin  A 

[  2  r  2  2  4 v  r  '  7 

1—  —  3  _  _  3__  1 

~-{3Ur<ll  +  Utq2)cosA  +  -(Urq2+Utq1)sm3A+-(Urql-Uiq2)cos3Ay-Q.ql  cos  i 

PkC  —  —  —  —  —  — 

A  =  n  + - B-{-'i(u.q.  +Uq  )  +  4£7  sin  A +  4(7  cosA  +  ((7  a  -Uq2  )cos2A 

2«a  1  t  /  t  r  \  t  ) 


(Urq2  +  £/g])sin2A}-Qcosi 


+ 


c*  =  o 


(109) 


Averaging  Eq.  (109)  by  the  mean  argument  of  latitude  over  one  orbit  and  ignoring  the  shadow 
effects  yields  the  long  period  and  secular  rates  [28], 
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-  3PkC„  - 

i  = - S-U.q, 

2  na  hl 

A  3 PkCR  - 

Q  = - —U.q 

2  na  sin/ 

3 PkCn  - 


3  PkCB  v 


-  3  PkCR  f-  -  , 

cB  =  o 


(110) 


3  cos/ 

2  na  sin/ 


PkC„U,a2 


32.2.2.2  Short  Period  Terms 

According  to  the  method  of  averaging  the  zero  order  terms  of  the  right  sides  of  Eq.  (109)  should 
keep  the  same  formulations,  but  in  terms  of  the  mean  elements.  Subtracting  Eq.  (110)  from  the 
zero  order  expansion  of  Eq.  (109)  and  integrating  the  remaining  terms,  the  short  period  terms  are 
[28] 


2  pur  r  —  —  i 

8a  = - £/..  cos2  +  t/.sin2  +  - 


Ur  cos  2  +  Ut  sin  2  +  {Urql  -  Utq2 )  cos  2A  +  -^(Urq2  +  Utqx )  sin  2A 


2  n2  a 


-  2  sin  A  +  -^qt  sin  2 A  ~~q2  cos  2 A 

-  2  cos  A  +  —  qx  cos  2 A  +  —  q2  sin  2 A 
/ 1_  2  2 


PkC  —  —  —  —  —  — 

8qx  = - jZ-  { Ut  sin  2 A  -  Ur  cos  2 A  - ( Urqx  +  3 U,q2 )  cos  A  +  ( 5 Urq2  -  3Utqx )  sin  A 

4  n  a  K  ' 

~(~Urqt  +  Utq2  j  cos  3 A  +  [Urq2  +  Utqx  jsin3/l|  +q2  cos  idPl 

PkC1  —  —  —  —  —  — 

Sq2  = - cos  22  +  Ur  sin  22  -  ( 3 Urq2  -  5 Utqx )  cos  2  -  ( 3Urq[  +Utq2 )  sin  2 

4  n  ay  v  ’  v  ’ 

~{Urq2  +Ulql^cos3A  +  (Urql  -f/(r/2)sin32|  —  qx  cos idPl 

T)J 

8A  =  — -T^-\Ur  sin  A-Ut  cos  2  +  (Urqx  -Er^2)sin22-(t7f.g2  +£/,g1)cos22]-$T2cos/ 
n a  L  v  /  v  /  j 

8Cb=  0 


(111) 


The  elements  on  the  right  sides  of  Eqs.  (110-111)  are  the  averaged  or  mean  elements. 
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3. 2.2.3  Extended  GA  STM  Including  SRP  Perturbations 

3.2.2. 3.1  Differential  Mean  Elements 

The  STM  or  the  partial  derivative  matrix  of  the  current  orbit  elements  with  respect  to  the  initial 

orbital  elements  in  the  mean  space  is  represented  by  - - •  Assume  the  secular  and 

0eo 

long  period  rates  are  constant.  The  mean  propagation  equations  are 

a  =  aQ  +  d(t  —  tQ ) 

A  =  A0  +  X(t-t0) 

i  =  i0  +  J(t-t0) 

=  01O +  ?!  (*“*<> )  (112) 
q2  =  q20+q2{t-to) 

a  =  a0+a(t-t0) 
cB  =  cB0  +  cB(t-t0) 

The  expressions  for  (pe  are  presented  in  Appendix  C. 

3.2.2. 3.2  Mean  to  Osculating  Transformation 

The  osculating  elements  are  obtained  as  functions  of  the  mean  elements 

e  =  e  +  Se  (113) 

where  e  are  the  mean  elements  and  Se  are  the  short  period  terms  as  shown  in  Eq.  (112)  and  Eq. 
(113),  respectively.  Taking  derivatives  with  respect  to  the  mean  elements  for  Jacobian  matrix,  we 
have  the  mean  to  osculating  transformation  matrix  or  D  matrix,  as  shown  in  the  Appendix  C. 

3.2.2. 3. 3  Geometric  Transformation 


The  Sigma  matrix,  X,  is  used  to  transform  osculating  elements  to  LVLH  coordinates.  Although 
the  solar  radiation  pressure  or  atmospheric  drag  is  a  nonconservative  perturbation,  the  velocity 
components  have  the  same  formulations  as  those  of  two-body  problem.  The  normal  accelerations 
are  given  in  Eq.  (105). 

3.2.3  Dynamic  Model  Expansion:  Drag  Perturbations 

Atmospheric  drag  is  one  of  the  major  perturbations  that  influence  satellite  motion,  especially  for 
satellites  in  low  earth  orbit  (LEO).  Extensive  modeling  has  been  accurately  built  up  for  the 
gravitational  field,  but  analytical  treatment  of  atmospheric  drag  is  very  difficult  since 
atmospheric  density  is  fluctuated  and  tabular  due  to  solar  activity.  It  is  well  known  the  dynamic 
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system  with  drag  is  nonlinear  and  nonconservative  and  the  analytical  solutions  of  the  dynamic 
system  are  unavailable.  The  best  we  can  do  is  to  derive  approximate  analytical  solutions  with 
desired  accuracy. 

Brouwer  and  Hori  [29]  obtained  the  analytical  solutions  for  gravity  and  drag  perturbations  using 
von  Zeipel  transformations  in  Brouwer  [3]  on  the  basis  of  an  exponential  density  model.  The 
density  function  was  expanded  into  a  series  of  the  ratio  of  the  eccentricity  to  the  density  scale 
height.  The  resulting  theory  is  limited  to  Earth  satellite  orbits  with  small  eccentricities.  Lane,  et 
al.  [30]  improved  the  Brouwer-Hori  work.  By  using  a  power  law  density  model,  the  expansion  of 
the  density  function  was  avoided  so  that  a  better  convergence  for  low  perigee  heights  is 
achieved.  Hoots  [31]  further  improved  Lane's  theory  by  eliminating  small  eccentricity  divisors  in 
the  differential  equations.  Moreover  he  applied  an  averaging  method  so  that  the  differential 
equations,  rid  of  the  mean  anomaly,  could  be  integrated  more  easily. 

The  differential  drag  has  been  used  for  controlling  relative  positions  between  satellites  [32-33], 
For  formation  flying,  the  difference  between  the  effective  cross  sectional  areas  leads  to  the 
differential  drag,  either  by  change  satellite  attitude  or  by  using  deployable  drag  panels.  Most 
papers  about  differential  drag  for  formation  control  are  based  on  CW  equations  or  more 
accurately  based  on  Schweighart  and  Sedwick  equations  [34]  in  which  the  differential  drag  is 
projected  into  the  LYLH  frame  or  along- track  direction. 

Reid  and  Misra  [35]  revised  the  Schweighart  and  Sedwick  equations  to  include  differential  drag 
and  extend  the  equations  for  reference  orbits  with  small  eccentricity.  Recently  Kumar,  et  al  [36] 
examined  the  feasibility  of  formation  maintenance  using  environmental  forces,  especially  solar 
radiation  pressure  and  aerodynamic  forces.  It  is  assumed  that  the  satellites  are  equipped  with 
solar  flaps  or  aerodynamic  flaps.  By  appropriate  rotation  of  these  flaps,  it  is  possible  to  influence 
the  relative  motion  between  satellites  in  a  formation. 

Most  authors  considered  the  differential  drag  as  a  means  to  control  satellite  formations.  Mishne 
[37]  introduced  mean  rates  due  to  drag  into  relative  motion,  which  is  a  basis  for  control.  Carter 
and  Humi  tried  to  include  the  aerodynamic  forces  into  the  analytical  solutions  of  relative  motion. 
Carter  and  Humi  modified  the  Clohessy- Wiltshire  equations  to  include  the  effects  of  atmospheric 
drag  in  two  separate  papers.  The  relative  motion  equations  developed  included  a  drag  force  that 
was  proportional  to  linear  velocity  [38],  Then  the  relative  motion  equations  developed  included 
the  effects  of  a  more  realistic  drag  model,  one  where  drag  is  proportional  to  the  square  of 
velocity  [39],  Based  on  these  simplifying  assumptions,  a  set  of  linear  differential  equations  were 
obtained  which  can  be  solved  in  terms  of  integrals.  This  enables  the  representation  of  the 
solution  of  the  problem  in  terms  of  a  state-transition  matrix.  Palmerini  et  al.  [40-41]  showed  the 
performance  based  on  the  Carter  and  Humi  solutions  is  more  efficient  than  CW  equations  if  the 
aerodynamic  perturbations  are  modeled  into  the  analytical  solutions. 

In  this  report  we  extend  the  GA  STM  [1]  to  include  differential  drag.  The  exponential 
formulation  is  selected  as  the  density  model.  The  method  of  averaging  is  used  for  setting  up  the 
STM  for  mean  element  propagation.  The  short  period  terms  are  ignored.  The  extended  GA  STM 
includes  parameter  variables  that  allow  modeling  different  ballistic  coefficients. 
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3. 2. 3. 1  Drag  Perturbations 

The  satellite  accelerations  due  to  atmospheric  drag  are  given  by 

XJ  =  -\CDCBpv\  (114) 

Assume  a  non-rotating  atmosphere,  the  acceleration  in  the  normal  direction  is  zero  and  the  radial 
and  tangential  components  are 

Ur  =  -\kpvvr 

U,  =  -\kpvvt  (115) 

£4=0 

where 


v,-=^esinf 
v<  =\[i(l  +  ecosf) 
v  =  'Ji{l+2ecosf+e2)/2 


and  p  is  the  atmospheric  density,  v  is  the  speed  of  the  satellite  relative  to  the  atmosphere,  v  is  a 
velocity  vector  in  the  direction  motion,  Cd  is  the  dimensionless  drag  coefficient. 

The  atmospheric  density  model  in  Eq.  (114)  is 


P  =  Ppe 


Ar~rp) 


(116) 


where  rp  is  perigee  altitude  and  /?  is  the  inverse  of  the  atmospheric  scale  height. 

3.2.3.2  Secular,  Long  and  Short  Period  Terms  from  Drag  Perturbations 

Since  the  normal  perturbed  acceleration  is  zero,  the  extended  nonsingular  Gauss  equation 
becomes 
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da 

dt 

d) l 
dt 

di 

dt 

jgi 

dt 

d<h 

dt 

dQ 
dt 
dC B_ 
dt 


2a2 


(qt  sin 6-  q2  cos 0)U r  +  —  1/ 


=  ft-- 


£/ 


J_ecos/ 

V  a  \+q 


ri  (  r  ) 

l+- 

V1  l+T/'v  /V 


esin  / 


=  0 


-hin«/  + 

"I  ' 


v  pj 


cos  0  -\ —  q. 
P 


U. 


=  —<-cos0U  + 


=  0 


Yf  \ 

\+r- 

r 

sin  0  +  —  q. 

u] 

\  p) 

P 

J 

=  0 


(117) 


Assuming  small  eccentricity,  we  have  approximately 

6  =  /l  +  2qx  sin  2  -  2q2  cos  2 
~  =  ^~ch  cos  X-q2  sin  2 

a  (118) 

sin  0  =  -q2  +  sin  2  +  qx  sin  2/1  -  q2  cos  2/1 
cos  6  =  -qx  +  cos  /l  +  qx  cos  2 2  +  q2  sin  22 

With  small  eccentricity  assumption,  expanding  Eq.  (1 16)  to  the  second  order  by  Taylor  series 


P  =  Pp  exp(-^e) 


\  +—%2e2  +  %{qx  cos  X  +  q2  sin  2)  +  —  p2  (ql  -q22  )cos22  +  —  %2q1q2  sin  22 


(119) 


where  d,  -  Pa  . 

Substituting  Eq.  (115)  and  Eqs.  (118-119)  into  (117)  gives 
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a  =  -na2CDCBpp  exp(-<i;e)[l  +  ^-£2e2  +(£,  +  3)(^  cos  h  +  q2  sinA)  + 

) cos  2?l  +  ^  sin  2^\ 

X  =  n  +  ^naC^C^p^exp^-^e)^  sin  k-q2  cosA) 


i  =  0 


qx  -  --naC DC Bp  pe*v{-&)[{\  +  £)qx+  2  + -E,2  q\  + -E,2  q\ 
^■i %2qlq2  sin  A  +  (E,  +  3)^  cos2A  +  q2  sin2A)  + 

^  <f  [q\  ~  q\ ) cos 3^  +  ^ Vi^2  sin 3^] 

=  --naCDCBp/)exp(-^)[(l  +  ^)g2+  2  +  -<fg2  +  -<fg2 

^ 2qlq1  cos  A  +  +  3)(^j  sin2A  -  q2  cos2A)  + 

7  <f  ~  ^ ) sin 3^  -  \  Z 2 ^2  cos 3^1 


cos  Ad- 


sin  A  + 


4 

Q  =  0 

c,  =  o 


(120) 


Averaging  Eq.  (120)  over  the  mean  argument  of  latitude  over  one  orbit,  we  obtain  the  long 
period  and  secular  rates 


a  =  -na2CDCBpp  exp(-f e)f  1  +  -t2  e 


a  3  ti 

A  -  n - a  At 

2  a 

i  =  0 


q1=-^naCDCBppexp(-%e)(l  +  %)q1 
q2=-]-nacDCBppexp(-&)(l  +  £)q2 


(121) 


£2  =  0 

c*  =  o 
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The  zero  order  terms  of  the  right  sides  of  Eq.  (120)  should  keep  the  same  formulations,  but  in 
terms  of  the  mean  elements.  Subtracting  Eq.  (121)  from  the  zero  order  expansion  of  Eq.  (120) 
and  integrating  the  remaining  terms  the  short  period  terms  are 


8a  =  -a2CDCBpp  exp(-^e)[(^  +  3)(^!  sin2 -q2  cos/l)  + 

£  [ql  ~  q\ )  sin  2A~  ?qxq2  cos  2/1] 

1  3/2 

8/ l  =—aCDCBpp  exp(-<^e)(^1  cos/l  +  q2  sin2)-  —  J— 8adt 


2 

8t  =  0 


1 


1 


Sqi=--aCDCBppexp(-te)[  \2  +  -fqf+-fq22 


sin2- 


^ 2qxq2  cos2  +  -^-(£ +  3)^  sin 22.  —  q2  cos 22.)  + 
T  f 2  (it  ~ ll ) sin 3/1  - ^ f29i92  cos 3/] 
S%=-\aCDC,pt  exp(-^)[-f2  +  ^-  +| 

sin2--^-(<f; +  3)(<^  cos 2A  +  q2  sin  22)- 


cos2  + 


12 

<5Q  =  0 


—  <^2  (q2  -ql )cos32  -  —  %2qlq2  sin32] 


8Cb  =  0 


(122) 


The  long  period  and  secular  rates  Eq.  (121)  are  obtained  through  the  second  order  expansion  of 
the  density  function;  however,  a  more  accurate  method  to  obtain  the  mean  rates  is  to  use 
modified  Bessel  functions  of  the  first  kind  [14]: 


a  =  —a2nkppf1 
e  =  -ankrfppf2 


i  =  0 
Q  =  0 
£0  =  0 

M  =  nn  — 


1% 
2  a„ 


-at 


(123) 
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where 


f,  =  [/„  +  2  el,  +  ;«“(/„  +  2. )]  exp  (-fiae) 

A  =  [2,  +  Hh  +  C)+W  (H,  +  A)]exp(-/?ae) 

In(/3ae)  are  modified  Bessel  functions  of  the  first  kind, 


or 


7«  (z) = 


OO  1 

y _ 1 _ 

f^(n  +  k)\k\ 


\n+2k 

<2; 


(124) 


(125) 


3. 2.3. 3  Extended  GA  STM  Including  Drag  Perturbations 
3.23.3. 1  Extended  GA  STM  for  Differential  Mean  Elements 

It  can  be  shown  that  the  solution  difference  is  very  small  between  Eq.  (121)  and  Eq.  (123)  for 
small  eccentricity.  Here  we  use  Eq.  (123)  to  set  up  the  mean  STM.  Assume  the  constant  mean 
rates 


a  =  a0  +  aDAt 


A  —  Aq  + 


f  3  n 

n - aDAt 

v  4  a  D 


At 


l  =  la 


qx=e  cos  co  =  (e0  +  eAt )  cos  a>0  = 

q2  =  e  sin  co  =  (e0  +  eAt )  sin  co0  = 
Q  =  Qn 


6  e  A 

1  +  — A  t 

V  eo  J 

W 

V  e0  J 


<ho  =(l  +  /A/)g,o 

#20  =  (l  +  AAt)^20 


(126) 


c  =  c 


where 
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(127) 


/3  =  -anCDCBT)2p  £- 


From  Eq.  (125),  we  have 


(128) 

We  can  remove  the  zero  eccentricity  in  the  function /? 

/3  =-^a«C£)Cfi^2/?p[(l+a^)/0 +(l-a^)/2 +|e(3/j +/3)]exp(-^ae)  (129) 

Taking  the  derivative  of  Eq.  (126)  with  respect  to  the  initial  mean  elements,  we  have  the  mean 
STM,  as  shown  in  Appendix  D. 

3 .2.3 .3 .2  Mean  to  Osculating  Element  Transformation 

Taking  the  derivative  of  Eq.  (122)  with  respect  to  the  mean  elements  yields  the  Dd  matrix,  as 
shown  in  Appendix  D. 


3.2. 3. 3. 3  Geometric  Transformation  Matrix 

Since  Uh  ,  the  normal  acceleration  caused  by  drag  perturbations,  is  zero,  substituting  in  Eq. 
(126)  gives  the  geometric  transformation  matrix. 

3.23.4  Drag  and  Gravity  Combined  Perturbations 
3.2. 3.4.1  Extended  GA  STM  for  Mean  Elements 

The  combined  mean  element  propagation 


a  —  aQ  +  aDAt 


A-T0  + 


'  .  3n  ^ 

n  +  ?lG~^  ~aDAt 


At 


<h  =  (! ' +  AAt)(<ho  cos(  A©g  )  -  q20  sin(  AcoG )) 
42  =  (! ' +  f3*d){q  io  sin(  Amc )  +  q20  cos(  A  coG )) 
Q  =  Q  +  Q 

u  o 

c  =  c 

^ B  ^B0 


(130) 


Taking  the  derivatives  of  Eq.  (130)  with  respect  to  the  initial  mean  elements,  gives  the  mean 
STM  for  drag  and  gravity  combined  perturbations. 
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3.23.4.2  Mean  to  Osculating  Element  Transformation 
The  mean  to  osculating  transformation  is 


D  =  Dg  +  Dd 


(131) 


The  drag  perturbations  do  not  change  the  geometric  transformation  matrix. 


3. 2. 3. 5  Correction  for  Coupling  Effects  from  Combined  Perturbations 

The  complexity  of  the  oblateness  and  drag  interactions  was  illustrated  by  Brouwer  and  Hori  [30], 
Green  [41]  expanded  the  combined  perturbations  to  the  second  order  and  considered  the  crossing 
terms  to  be  the  coupling  effects.  The  simple  way  to  deal  this  problem  is  to  use  osculating  height 
for  air  density  when  calculating  mean  drag  functions  [10,  42-43].  This  is  because  the  osculating 
height  variations  due  to  the  Earth  oblateness  for  LEO  satellites  can  be  several  kilometers,  which 
greatly  influence  the  atmospheric  density. 


The  correction  for  the  height  or  radial  distance  [4]  is 


Re 

( 

3  .  2 

\ 

r 

ecos  f^ 

1 

— sin 

i 

1  +  2  —  + 

p 

2 

) 

art 

1  +  9  ) 

— J,  — -  sin2  i  cos  29 
4  2  p 


(132) 


Notice  the  variables  on  the  right  sides  of  Eq.  (132)  are  the  mean  variables.  The  corrected  density 
is 


P  =  PPe 


P(r+Sr~rp ) 


(133) 


Substituting  Eq.  (132)  into  (133)  then  Eq.  (130)  gives  the  mean  rates  for  the  combined 
perturbations. 

3.2.4  Semi-Analytic  Extended  GA  STM  Including  Drag  and  Gravity  Perturbations 

Accurate  modeling  of  drag  perturbations  can  be  achieved  by  the  use  of  semi-analytic  schemes. 
Arsenault  et  al  [44]  suggested  using  a  lower-order  numerical  integration  scheme  for  drag 
perturbations.  Hoots  [10]  applied  the  Gauss-Legendre  formula  to  integrate  the  averaged  drag 
functions. 


Since  the  short  period  terms  are  averaged  out,  we  use  the  lower  order  integration  methods  to 
propagate  the  mean  elements.  The  Euler  method  is  an  example: 

xn=xn-l+f(xn-l’tn-l)h  (134) 
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where  h  is  the  step  size.  The  GA  STM  semi-analytic  propagation  is  achieved  using  the  steps 
outlined  by  Eqs.  (135-137).  The  differential  orbital  elements  at  the  initial  time  t()  are  obtained 
from  the  initial  relative  state  vector  as 

^e(f0)  =  Zr1(f0)E-1(f0)X(f0)  (135) 

The  averaged  elements  of  the  chief  and  deputy  are  propagated  for  an  arbitrary  time  period 
[t0 ,  tN  ]  by  using  the  Euler  method  with  n  steps 


eC  (4  )  =  eC  (4-1  )  +  /'(ec  (4-1  )’4-l 
eD  (4  )  =  eD  (/«- 1  )  +  f  (eD  (4-1  )  » 4-1  )  ^ 

ec  (4)  =  e(4 ) 

(4) =  e(4)  +  <^e(4) 

ec  (^V  )  =  eC  (4  ) 

eD  (4f)  =  eD  (4) 

Finally,  the  relative  state  vector  at  the  desired  final  time  is  obtained  as 

^■(4r)  —  ^-'(4r)^(4r)(eo  (6v)_ec  ((v)) 


(136) 


(137) 


where  e  is  the  vector  of  the  averaged  elements  and  the  subscripts  C  and  D  indicate  the  chief  and 
deputy  satellites,  respectively.  An  alternative  option  is  to  use  integration  methods  with  variable 
step  sizes,  such  as  ODE45  in  MATLAB.  The  averaged  semi -major  axis  and  eccentricity  values  are 
obtained  by  the  Legendre-Gauss-Lobatto  integration  rule  over  one  orbit  using  the  rates 


ad=-^^Pi(1  +  ecos  Ei  )K  e  cos  Ei  Y1  wi 

2  l=\  (138) 

ed=~  Z  Pi  0  +  e cos  Et  )K  ( 1  - ecos E .  \Yl  cos E,wt 

1  i= 0 


p  is  the  atmospheric  density  model,  which  in  general  is  represented  by  a  complex  model,  such  as 
a  tabular  data  model,  here  the  exponential  model  is  used  as  an  example. 

3.3  Dynamic  Model,  Navigation  Accuracy  and  Thruster  Accuracy  Consistency 

Although  the  guidance,  navigation,  and  control  (GNC)  subsystems  are  commonly  grouped 
together  in  the  spacecraft  design  process,  in  practice,  there  are  still  inefficiencies  caused  by  a 
lack  of  interdependency  between  them.  Notably,  the  dynamical  model  in  the  navigation 
algorithm  is  often  set  ad  hoc  without  explicit  regard  for  the  level  of  measurement,  guidance,  or 
control  errors  expected.  We  understand  qualitatively  that  the  “best”  dynamical  model  will  meet 
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some  user-defined  criteria  on  state  uncertainty  and  maneuver  cost  while  minimizing 
computational  effort.  For  example,  it  is  operationally  suboptimal  to  implement  a  high-fidelity 
gravity  model  when  the  measurement  hardware  is  known  to  be  wildly  inaccurate  [46].  On  the 
other  hand,  maneuvers  executed  to  cancel  out  the  effects  of  known  yet  ignored  perturbations,  if 
frequent  enough,  would  deplete  on-board  fuel  very  quickly  and  shorten  mission  lifespan.  A 
quantitative  study  is  thus  warranted,  preferably  via  an  analytical  approach,  on  the  consistency  of 
GNC  hardware  subsystems  versus  dynamical  modeling  algorithms. 

In  this  chapter,  we  tackle  a  subset  of  this  problem:  namely,  we  develop  methods  to  quickly 
survey  the  trade  space  between  navigation  system  parameters  and  dynamical  model  fidelity.  We 
focus  our  efforts  on  forces  that  have  precise  deterministic  physical  models,  e.g.,  the  Earth's 
gravity,  such  that  modeling  errors  may  be  regarded  as  biases.  Our  contributions  are  as  follows. 
First,  we  show  that  a  change  in  the  dynamical  model  may  be  related  to  its  corresponding  change 
in  STM  via  an  explicitly  defined  vector  function.  This  result  not  only  allows  us  to  compute  the 
state  transition  matrix  for  multiple  dynamical  models  efficiently,  but  also  demonstrates  that  the 
change  in  the  STM  is  typically  on  the  order  of  10°  %  over  tens  of  orbit  periods  of  the  chief 
satellite.  As  such,  as  our  second  contribution,  the  linear  sensitivity  relating  dynamical  model 
fidelity  to  maximum  a  posteriori  state  estimate  bias  is  derived.  Although  the  dynamical  model 
error  cannot  be  factored  out  as  a  linear  operator  like  the  observation  covariance  or  cadence,  we 
may  still  gain  approximate  yet  quick  design  insight  into  choosing  an  appropriate  dynamical 
model  for  a  set  of  given  navigation  system  requirements.  Finally,  the  cost  of  employing  a 
particular  dynamical  model  will  be  characterized  through  a  maneuver  metric  first  proposed  by 
Schaub  and  Alfriend  [47],  That  is,  the  state  deviation  of  the  estimated  trajectory  from  the 
reference  is  analytically  transformed  into  a  Ay  via  Gauss'  planetary  equations.  The  three  results 
presented  will  simplify  the  workflow  of  designing  GNC  systems  by  mitigating  the  need  to 
conduct  a  large-scale  numerical  validation  of  system  performance. 

3.3.1  Relating  Dynamical  Model  and  State  Transition  Matrix  Fidelity 

Suppose  x{t ;  x0)  is  the  trajectory  given  by  a  high-fidelity  dynamical  model  x(t;  x0)  that  we 
regard  as  truth,  and  xF(t;  x0)  is  that  given  by  the  approximated  dynamics  xF(t;  x0).  The 
semicolon  separates  the  time  variable  t  from  the  initial  conditions  x0,  which  are  treated  as 
parameters.  We  define  a  scaling  function  m(t)  such  that 

m‘  (0  xj,(f;  xo)  =  x*'(t;  x0)  (139) 

for  all  t  and  all  i,  where  the  superscript  indicates  a  vector  or  matrix  index.  We  assume  that  the 
vectors  are  expressed  in  an  inertial  reference  frame. 

The  state  transition  matrix  (STM)  of  the  true  dynamics  is  given  by 

4>«((;x„)  =  A  (140) 

K 


Then,  for  O'T  ^  0, 
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(141) 


O'-'  dx'  /  dxJ0 

Op  dxydx^ 


x'(f;  x0  +  eO  -  xf(l;  x0)]  /s  _  x'(l;  x0  +  e7)  -  xf(l;  x0) 
xj.(f;  x0  +  e^)  -  x^(f;  x0)]  /e  xp(^  xo  +  e-0  -  x0) 


where  e7  is  a  small  variation  along  the  /-direction  with  magnitude  e«  1.  Now, 

x'(i;  xo)  =  f  x'(t;  x0)dr  +  x[,  =  f  m'(r)x|,(T,  x0)dr  +  xj,  (142) 

Jo  Jo 


Applying  integration  by  parts, 


:'(r;  x0)  =  [ni'(r)x^(r;  x0)£  -  m'(T)x^(r;  x0)dr  +  x‘0 


(143) 


tti1(t)xf(t,  x0)dr  «  0  is  a  good  assumption  to  order  10°  %  or  better  over  10°  days  for  many 
Earth  orbiters.  Therefore, 


x'(i;  X0)  «  m''(0xF(t;  x0)  -  m'(0)xj)  +  xl0 


(144) 


m( 0;  ;c0)  ~  m( 0;  x0  +  eJ)  «  1  is,  again,  a  good  assumption  to  order  10°%  or  better.  Thus,  we 
can  rewrite  Eq.  (141)  as 


G>iJ 


m'(f)  [x^(l;  x0  +  ej)  -  x‘F(t;  x0) 
xp(t;  xo  +  e-0  -  xj,(i;  x0) 


=  m'(/) 


(145) 


As  such,  one  may  compute  the  STM  of  the  approximated  dynamics  Of(1;  xo)  solely  through  the 
solution  flow  of  the  true  dynamics  x(l;  xo ),  its  corresponding  STM  0(1;  xo),  and  the  function 
m(l),  which  is  readily  available  via  xf(1;  xo).  That  is,  for  some  1  such  that  mfl)  ^  0, 


Oi7(Tx0) 

m'(t) 


(146) 


Since  there  is  no  need  to  solve  ODEs  for  each  component  in  Of,  the  method  above  significantly 
speeds  up  the  propagation  of  the  secondary,  assuming  it  is  within  the  linear  dynamical  regime  of 
the  primary.  Furthermore,  the  scaling  function  m(l)  provides  insight  into  the  accuracy  of  the 
approximated  STM  with  respect  to  the  truth.  The  difference  between  the  true  and  approximated 
STMs  AO  is  given  as 


A 0(1;  xo)  =  0F(l;  x0)  -  0(1;  x0)  =  [M(l)  -  I]  0(1;  x0) 
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(147) 


where  M{t )  is  a  matrix  with  m (t)  along  its  diagonal  and  I  is  the  identity  matrix. 
For  instance,  consider  the  simple  1 -dimensional  case 


x  =  x 


xF  =  mx 


(148) 


where  m  is  held  constant.  Then,  the  solution  flow  for  the  true  dynamics  is  given  as  x(t;  x0)  = 
x0et_t°.  Next, 


xF  =  I  mx0eT  f°  dr  +  x0  =  mx0et  to  =  mx 
to 


(149) 


such  that  <t>F  =  mO  as  derived. 

3.3.2  Linear  Sensitivity  Analysis  Between  Dynamical  Model  Fidelity  and  Maximum  A 
Posteriori  State  Estimates 

We  now  derive  below  such  a  function  for  the  batch  least  squares. 

Nominally  [48], 


x0  =  A-,N 
Po  =  A"1 

N 

A  =  +  P0 

1 

N 

N  =  2(//7P71y;)  +  p-1x0 

i 

Hi  =  HMttx o) 


(150) 


(151) 


(152) 


where  x0  and  Po  are  the  state  deviation  estimate  and  associated  covariance  at  epoch,  respectively, 
A  is  the  information  matrix,  N  is  the  normal  vector,  Pi  is  the  observation  covariance,  P0  is  the  a 
priori  covariance,  y/  is  the  observation  deviation  vector  from  the  reference  trajectory,  Ht  is  the 
linear  observation-state  relationship,  and  0(tj;  x0)  =  <J>j  is  the  STM.  Subscripts  indicate  values 
for  the  z'-th  observation,  and  the  summation  is  over  all  zV  observations.  We  assume  that  the 
observations  are  uncorrelated  in  time;  thus,  both  the  information  matrix  and  normal  vector  may 
be  accumulated  as  a  summation. 

Now,  suppose  we  have  an  approximate  dynamical  model  such  that  the  new  trajectory  is  given  as 
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Then,  let 


xF(f/;x0)  =  x(ty;x0)  +  Ax(ti;xo),0F(to;xo)  =  O(t,-;x0)  +  AO(t;;x0) 
«  XF((  =  X;  +  AX;,  0Fii-  =  O;  +  AO; 


(153) 


xF  =  (Af)  1  Nf,Pf  =  (AF)  1  (154) 

denote  the  state  deviation  estimate  and  associated  covariance,  respectively,  obtained  from  the 
approximate  model.  The  goal  is  to  express  analytically  Ax0  =  xF  0  —  x0,  APo  =  Pf,o  -  Po.  First, 
simply  plug  in 

Ap  =  Yj  [(®i  +  AO ifHjRj'HiiOi  +  AO;)]  ~  A  +  £  [A HJRT1Hi  +  HJRt'AH,]  (155) 


where  A Ht  =  Ht  AOj  and  any  second-order  terms  are  ignored.  Similarly, 

Np  =  £  +  AOifHTp-Vr.i]  +  Po^o 


Now,  given  yFi  =  yt  +  A  yt  ~yi~  Ht  Axt  since  Yt  -  G(xFi)  =  Yt-  G(xt )  +  A  yt  A  yt  = 

G(xt )  —  G(x  pj), 

Nf  =  N  +  £  [A HjRj'y;  -  HJRt'H.Ax,]  .  (157) 

So, 

(a  +  Yj  [a hJrj'Hi  +  H^'Atfi])  Xp,i  =  N  +  Yj  [a HjRj'yi  -  HjR~]H;Ax]  (i58) 


xF  =  [A  +  LYX  (N  + 1}) 


(159) 


where 


L  =  Yj  [a HjR-'Hi  +  HjRJlAHi] ,  //  =  £  [Afffa"1*  -  HjR^HAxi 


(160) 


Note  that  for  some  square  matrix  Q  such  that  Qn  ->  0  (n  — >  oo), 

(I  +  Qr'-I-Q  (161) 

Then,  for  some  square  matrices  A  and  B  with  the  same  dimensions  and  A  is  semi-positive 
definite, 
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(A+  B)_1  ~  A-1  -A_1BA_1 

if  (A^BA-1)"  — »  0  (n  —>  oo).  Substituting  A  =  A  and  B  =  L, 

xF  «  (A-1  -  A-,IA-1)(N  +  T})  »  PN  +  Pt}  -  PLPN 

where  P  =  A  1  and  again  ignoring  second-order  and  higher  terms.  So, 

(  Ax  =  Pi}  -  PLPN  =  Pi}  +  APN 
1  A  P  =  ~PLP 


(162) 


(163) 


(164) 


These  expressions  show  that,  unlike  the  measurement  covariance  matrix  (if  assumed  constant 
over  all  observations)  or  the  measurement  cadence,  the  effects  of  STM  bias  cannot  be  factored 
out  from  the  normal  equation.  Nonetheless,  if  we  substitute  Eq.  (159)  into  the  expression  for  Af, 

=  Z  [tfMiHjRi'HiMi <Di]  (165) 


But  1HiMi  <  (mf)2Hj  Rt  1Hi  where  mf  is  the  component  of  m(t,)  =  rtij  with  the 

maximum  absolute  value.  Thus, 


Af  —  A  < 


mM  ~  1  ]  [mM  + 


(166) 


where  m%  is  the  maximum  of  mf  over  all  observations.  If  —  1  is  on  the  order  of  10'3  as  in 
Fig.  43,  then  the  difference  in  the  information  matrix  between  the  two  models  is  also,  at  most,  on 
the  order  of  10"3.  In  conclusion,  a  small  bias  in  the  STM  has  an  equally  small  effect  on  the 
information  matrix,  and  thus  the  estimated  covariance. 


Although  we  omit  the  details  of  the  derivation  for  brevity,  we  may  derive  similar  results  for  the 
Kalman  filter  such  that,  for  the  z'-th  observation, 

Ax,*  =  Ax,*  +  AKiiyi  -  Hxj)  -  K,H(Ax,  +  Ax,)  ^ 

A  Pi  =  (I  -  KiH)APi  -  AKjHPj  =  (I  -  K/Pf)AP,(/  -  KtH)T 


where 


Ax,-  =  A<M;_i  +  C>Ax;_i 
A  Pi  =  A  <hP,_,0T  +  OP/.,  AOt  +  (DA  P/_1Ot 
A Ki  =  (I  -  KiH)APlf{HPHT  +  R)~x 
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4  RESULTS  AND  DISCUSSION 


4.1  Gravitational  Perturbations 

4.1.1  Kinoshita  and  Kaula  Models 

Figure  1  shows  the  errors  in  the  three  components  as  well  as  the  r.m.s.  error  in  the  relative 
position  obtained  from  the  Kinoshita  theory  when  compared  with  a  numerical  simulation.  The 
effects  included  are  Ji,  augmented  with  those  due  to  the  mean  and  long  periodic  effects  of  J22,  J4 , 
and  J6. 


Figure  1.  Modeling  J2,  J22,  J4,  Jh  (mean  and  long-periodic  effects) 


Effects  of  the  J21  Short-periodic  terms 

Figure  2  shows  the  effects  of  including  the  short  periodic  contributions  from  J21.  It  shows  the 
errors  resulting  from  modeling  J2  (mean,  short  and  long  periodic)  and  including  J22  mean,  long 
and  short  periodic  effects.  Figure  2  shows  that  the  errors  are  uniformly  lower  in  all  the  three 
axes  with  respect  to  those  in  Figure  1.  Hence,  modeling  the  J22  short-periodic  terms,  albeit  for 
e=0,  is  beneficial  and  improves  the  accuracy  of  the  relative  motion  STM.  Figures  3-5  show  the 
differences  between  the  augmented  Kaula  and  GMAT  solutions  in  the  radial,  along-track  and 
cross-track  directions. 
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Figure  2.  Error  from  first  order  of  Ji ,  J4  and  Jt,  modeling  plus 
J2 2  mean  rate  and  short  period 
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Figure  3.  Difference  between  the  Kaula  and  GMAT  with  20x20 
gravity  field  without  J22  secular  and  short  period  terms 
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Figure  4.  Difference  between  the  Kaula  and  GMAT  with  20x20 
gravity  field  with  J2 2  secular  terms  and  without  J2 2  short 

period  terms 


x  10'4  x  10'3 


Figure  5.  Difference  between  Kaula  and  GMAT  with  20x20 
gravity  field  with  J2 2  secular  and  short  period  terms 
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Table  1.  Root  mean  square  error  for  the  Extended  GA-STM 
propagation  (10  days) 


Effects  modeled  for 
determining  the  initial 
conditions 

Modeling  Perturbations 
in  the  extended  GA  STM 
for  propagation 

ICs  obtained  from  the 

classical  Kaula  model 
(km) 

J2  only 

J2  only 

0.009433 

A  J22 

72,722 

0.006497 

J2 5  J2  5  J2  sp 

72,722 

0.004661 

*/2?  J2  ?  J2  Sp,  J3 

J2,  J22  and  J3 

0.006539 

J2,  722,  J22  sp,  7),  J4 

J 2?  J 2  3 1 J4 

0.005952 

t/2,  2  5  J2  Sp,  «/j,  J 4-,  J5 

J2  ?  J5 

0.006566 

7?,  J2 ,  7?  sp,  7?,  J 4t  7), 
J6 

«^2?  J2  ?  J 3i  J 4r>  J5i  J 6 

0.005540 

4.1.2  GA-STM  with  Hoots  Variables 

The  Hoots  version  of  the  GA-STM  can  be  evaluated  via  MATLAB  simulation.  Its  error  (versus 
a  numerically  integrated  trajectory)  can  be  compared  with  the  results  from  directly  differencing 
the  analytically  propagated  trajectories  of  both  chief  and  deputy  satellites.  These  numerically 
computed  differences  can  then  be  mapped  into  relative  space,  forming  a  rough  lower  bound  on 
the  error  for  closed-form  relative  motion  approximations  (such  as  the  GA-STM). 

In  Figure  6,  the  satellites  are  in  a  near-PCO  (Projected  Circular  Orbit)  formation:  the  zero-drift 
constraint  on  8a  is  turned  off,  the  relative  orbit  size  is  1  kilometer,  and  the  deputy's  initial  phase 
angle  is  0.  The  chiefs  osculating  initial  conditions  are  semi-major  axis  7100  kilometers, 
eccentricity  0.01,  inclination  50  degrees,  and  all  other  classical  elements  0.  J2  is  considered  as 
the  only  perturbation.  The  relative  position  and  velocity  error  magnitudes  are  shown  for  the 
Hoots-element  GA-STM  as  well  as  for  direct  differencing  of  first-order  orbit  theories  due  to 
Hoots  [10]  and  to  Kinoshita  [6].  In  this  scenario,  the  STM  shows  error  characteristics 
competitive  with  the  direct  differencing  models,  especially  for  position  error. 
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_ number  of  orbits _ 

Hoots  - Kinoshita  Hoots  GA  STM 


Figure  6.  First-order  theory  comparison:  Direct  differencing 
vs.  GA-STM 


The  new  version  can  also  be  compared  against  previous  versions  of  the  GA-STM  formulated  in 
terms  of  nonsingular  elements  and  equinoctial  elements.  Figure  7  shows  results  from  a  chief 
satellite  in  a  near-circular  low-Earth  orbit;  the  LVLH  components  and  the  magnitude  of  the 
relative  position  error  are  shown  for  all  three  STMs,  as  well  as  for  both  versions  of  direct 
differencing.  For  this  case,  the  STM  errors  are  not  significantly  greater  than  those  for  direct 
differencing. 
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Figure  7.  Relative  position  error  for  near-circular  orbit 


Figure  8  shows  the  same  relative  position  error  plots  for  a  different  scenario.  In  this  case,  the 
chiefs  eccentricity  is  approximately  0.1.  Here,  the  STMs  perform  slightly  worse  than  the  direct 
differencing  models,  not  an  unexpected  result. 
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Figure  8.  Relative  position  error  for  eccentric  orbit 


Figure  9  shows  the  relative  position  error  for  a  formation  in  geostationary  orbit.  Because  the 
chiefs  orbit  is  equatorial,  both  Kinoshita  theory  and  the  nonsingular-element  STM  are  singular 
and  are  not  shown. 
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Figure  9.  Relative  position  error  for  geostationary  orbit 


Note  that  in  many  cases  the  STM  errors  are  nonzero  even  at  the  initial  time.  In  order  to 
propagate  the  relative  orbit  accurately  with  an  STM,  the  relative  initial  conditions  must 
themselves  undergo  a  linearizing  transformation— this  transformation  introduces  an  offset  at  the 
initial  time,  but  prevents  an  accumulating  error  over  the  course  of  the  propagation. 

4.1.3  Numerical  Verification  using  GMAT  and  a  Graphical  User  Interface  (GUI) 

The  analytic  formulae  for  the  second-order  short-period  effects,  the  first-order  long-period 
effects  and  second-order  secular  effects  have  been  verified  for  the  zonals  J3  and  /4  with  the 
published  results.  A  MATLAB-based  code  for  the  relative  motion  STM,  including  all  the 
expressions  for  J2  effects  computed  by  using  Maple,  and  the  generalized  analytic  formulae  for 
the  higher  zonals,  has  been  developed.  A  graphical  user  interface  (GUI)  has  also  been 
implemented  to  easily  configure  and  run  simulations,  as  well  as  analyze  results.  An  interface  to 
National  Aeronautics  and  Space  Administration  (NASA)  GMAT  software  is  also  implemented  as 
part  of  the  GUI  in  order  to  compare  the  analytic  STM  accuracy  with  the  direct  numerical 
propagation  using  GMAT.  Figure  10  shows  the  screenshot  of  the  Matlab-based  Simulation  GUI. 
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For  comparing  the  accuracy  of  the  STM,  two  satellites:  chief  and  deputy,  are  simulated  using  the 
NASA  GMAT  software  with  70  x  0  JGM-3  gravity  model.  The  mean  initial  conditions  of  chief 
are  chosen  with  a=  7100  km  and  node  angle,  perigee  angle  as  well  as  the  mean  anomaly  all  as 
zero.  Two  different  reference  orbits  are  simulated:  eccentricity  0.01  with  inclination  50°  and 
eccentricity  0  with  inclination  0°.  The  mean  initial  conditions  of  deputy  are  computed  according 
to  a  projected  circular  orbit  type  formation  with  a  baseline  distance  of  1  km.  The  osculating 
initial  conditions  of  chief  and  deputy  required  for  numerical  propagation  are  computed 
analytically  using  the  mean  to  osculating  transformations  as  described  in  the  previous  sections. 
The  orbits  of  chief  and  deputy  are  propagated  for  ten  days  in  GMAT  and  the  relative  position  and 
velocity  states  in  the  curvilinear  frame  are  computed.  Using  the  extended  GA-STM,  the  same 
relative  states  are  also  directly  propagated  for  ten  days.  The  maximum  degree  up  to  which  zonal 
harmonics  are  included  in  the  extended  GA-STM  was  successively  increased.  Figure  1 1  shows 
the  absolute  position  and  velocity  root-sum-square  (RSS)  error  plots  for  chief  satellite  with 
eccentricity  0.01  and  inclination  50°  after  ten  days  of  propagation  as  the  maximum  degree  of  the 
extended  GA-STM  is  increased.  Figure  12  shows  the  relative  position  errors  for  a  PCO  type 
formation  with  1  km  of  baseline  distance  for  the  same  chief  reference  orbit.  It  is  noted  that  with 
the  maximum  degree  of  the  GA-STM  J20,  the  absolute  position  errors  were  reduced  to  less  than  1 
km  and  relative  position  error  less  than  1  m  after  ten  days  of  propagation.  Figure  13  and  14 
shows  the  absolute  and  relative  errors  when  the  reference  orbit  is  chosen  as  equatorial  and 
circular.  The  similar  trend  in  the  absolute  motion  as  well  as  relative  motion  was  observed  for  this 
case. 


Degree 

Figure  11.  Absolute  Position  and  Velocity  Error  vs.  Degree  of 
the  Extended  GA-STM. 
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Figure  12.  Relative  Position  Error  Versus  Degree  of  the 
Extended  GA-STM 
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Figure  13.  Absolute  Position  and  Velocity  Error  Versus  Degree 
of  the  Extended  GA-STM  for  Equatorial  Circular  Reference 

Orbit 
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Figure  14.  Relative  Position  Error  Versus  Degree  of  the 
Extended  GA-STM  for  the  Equatorial  Circular  Reference 

Orbit. 
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4.2  Non-Earth  Gravitational  Perturbations 
4.2.1  Third  Body  Perturbations 

A  representative  example  for  which  the  lunar  perturbation  is  significant  is  the  NASA 
Magnetospheric  Multiscale  (MMS)  mission,  designed  to  study  magnetic  reconnection,  charged 
particle  acceleration,  and  turbulence  in  key  boundary  regions  of  the  Earth’s  magnetosphere.  The 
MMS  orbit  is  highly  eccentric  and  lunar  perturbation  cannot  be  ignored.  The  starting  epoch  is 
August  5,  2013.  The  mean  inclination,  eccentricity,  and  right  ascension  of  the  lunar  orbit  are 
assumed  to  be  19.65°,  0.0497,  and  350.36°,  respectively.  The  lunar  position  data  at  30  minute 
intervals  is  taken  from  the  Jet  Propulsion  Laboratory  Horizons  website. 

The  reference  orbit  is  assumed  to  have  the  following  initial  mean  elements:  a=42095  km, 
e=0.81818, 1=28.5°,  D=357.857°,  00=298.2253°,  andMo=180°.  The  initial  phase  angle  ao  is 
selected  to  be  zero  and  the  relative  orbit  size  p,  is  10  km.  Because  of  the  high  eccentricity,  the 
relative  orbit  undergoes  an  expansion  near  perigee  and  therefore,  p  as  well  as  the  relative  orbit 
element  differences  do  not  remain  constant  along  the  orbit.  The  simulation  time  is  100  days. 
Earth's  gravitational  perturbations  are  not  included  in  the  simulations  presented  in  this  section. 
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4.2.1. 1  Results  for  the  Reference  Orbit 


The  equations  of  motion  for  the  chief  are  propagated  and  the  position  and  velocities  are 
converted  into  osculating  orbital  elements.  In  Fig.  15,  these  results  are  compared  with  the 
respective  solutions  obtained  from  Eqs.  (85). 
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Figure  15.  Averaged  Lunar  Model  vs.  Nonlinear  Model 


The  blue  lines  represent  the  averaged  model,  while  the  green  lines  stand  for  the  nonlinear  model 
in  Fig.  15.  One  can  see  that  there  are  obvious  differences  in  the  semi-major  axis  and  eccentricity 
for  the  two  models  due  to  the  use  of  the  same  initial  conditions  for  both  the  simulation  models. 
The  initial  conditions  for  the  nonlinear  simulation  should  be  corrected  to  fit  the  averaged  model. 
Application  of  the  least  squares  method  to  correct  the  initial  conditions  results  in  the  plots  in  Fig. 
16. 
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Figure  16.  Averaged  Lunar  Model  vs.  Nonlinear  Model  after 

Correction 


Figure  16  indicates  that  the  averaged  model  predicts  the  long-term  behavior  very  well,  providing 
confidence  that  the  effects  of  the  third-body  have  been  incorporated  correctly.  Only  the  first  term 
is  taken  in  the  infinite  series  represented  in  Eq.  (81)  for  the  averaged  potential.  From  Fig.  16,  we 
can  see  the  ( n=2 )  term  dominates  the  third-body  perturbation,  akin  to  the  gravitational  potential 
of  the  Earth,  in  which  Ji  plays  the  dominant  role.  Figures  15-16  also  show  that  it  is  important  to 
include  the  lunar  orbit’s  inclination  in  the  averaged  model  for  higher  accuracy,  as  concluded  in 
Reference  8. 

We  compare  the  nonlinear  simulation  results  with  GMAT  for  further  validation,  as  shown  in  Fig. 
17. 
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Figure  17.  Differences  between  Nonlinear  Simulation  and 

GMAT 


GMAT  is  an  open-source  space  mission  design  tool  to  model  and  optimize  spacecraft  trajectories 
in  flight  regimes  ranging  from  low  Earth  orbit  to  lunar,  libration  point,  and  deep  space  missions. 
GMAT  has  been  developed  by  a  team  of  NASA,  private  industry,  and  public  and  private 
contributors.  Note  that  only  the  lunar  perturbations  are  included  in  both  the  simulations. 

4.2. 1.2  Simulation  Results  for  Relative  Motion 

The  process  outlined  previously  for  the  simulation  of  the  reference  orbit  is  repeated  for  the 
deputy's  orbit.  In  this  section,  numerical  results  produced  by  the  extended  GA  STM  are 
presented.  To  isolate  the  third-body  effects,  we  neglect  the  Earth’s  gravitational  perturbations 
first.  Figures  18-19,  respectively,  show  the  errors  in  relative  position  with  and  without  modeling 
the  third-body  effect. 
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Figure  18.  Relative  Position  Errors,  Two-body  STM  vs.  Lunar- 
perturbed  Nonlinear  Model 


Figure  19.  Relative  Position  Errors,  STM  includes  Averaged 
Lunar  Perturbations 


Note  that  the  reference  orbit  for  this  example  has  a  large  semi-major  axis  and  is  also  highly 
eccentric.  Figure  18  shows  that  the  in-track  component  of  the  relative  position  error  vector  has 
the  largest  growth,  while  the  smallest  error  results  in  the  cross-track  direction.  Comparing  Figure 
18  with  Fig.  19,  it  is  seen  that  the  accuracy  greatly  improves  with  the  inclusion  of  the  lunar 
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perturbation  model  in  the  GA-STM,  especially  in  the  cross-track  direction.  The  residual  errors 
are  mostly  caused  by  absence  of  the  short  period  lunar  effects. 

The  Ji  perturbation  is  three  orders  of  magnitude  larger  than  that  due  to  the  higher-order  Earth 
gravitational  perturbations.  It  is  necessary  to  investigate  the  combined  effects  of  the  Ji  and  the 
third-body  perturbations.  In  the  examples  presented  for  error  comparisons,  the  initial  conditions 
for  the  numerical  integration  models  are  determined  by  the  least  squares  approach. 


Figure  20.  Relative  Position  Errors,  Two-body  STM  vs. 
Nonlinear  Simulation  With  Ji 


Figure  20  shows  the  position  errors  as  a  result  of  neglecting  Ji ,  arising  from  the  use  of  the  two- 
body  STM  model.  A  comparison  of  Figures  18  and  20  confirms  that  the  position  error  due  to  the 
Ji  perturbation  is  larger  than  that  caused  by  the  Moon;  although  the  semi-major  axis  of  the  MMS 
orbit  is  large,  about  42,000  km. 

Figures  21-22  present  the  improvements  obtained  from  the  use  of  the  extended  GA  STM  for 
propagating  relative  motion  and  for  the  least  squares  corrections  for  the  initial  conditions.  The 
extended  GA  STM  models  the  Ji  and  third-body  perturbations.  The  in-track  position  error  varies 
significantly  between  apogee  and  perigee,  proportional  to  the  size  of  the  relative  orbit. 
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Days 


Figure  21.  GA  STM  vs.  Nonlinear  Simulation  including  Ji  and 

Lunar  Perturbation 


Days 


Figure  22.  Extended  GA  STM  with  Lunar  Effect  vs.  Nonlinear 
Simulation  including  Ji  and  Lunar  Perturbation 


4.2.2  SRP  Numerical  Results 

The  mean  orbit  element  of  the  chief  satellite  is  a=43527. 7589899  km;  e=0.0005671;  Q=  18.88 
deg;  co=321.5388  deg;  M=38.4719.  The  cross  sectional  area  is  5.02  m2  and  the  satellite  mass  is 
611  kg.  The  simulation  time  is  10  days.  Only  the  SRP  perturbation  is  included  in  the  simulations. 
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The  inclination  i  is  set  as  1 1 3  deg  so  that  the  satellites  will  be  illuminated  all  the  time  to  avoid 
the  shadowing  effects. 

4.2.2. 1  SRP  Model  Verification 

Figures  23  and  24  show  the  accuracy  of  the  analytical  model,  with  and  without  corrections  for 
the  initial  conditions  of  the  analytical  model,  respectively.  The  solid  lines  stand  for  the  solutions 
from  Eq.  (107)  and  the  circles  represent  the  results  from  Eq.  (110).  Biases  in  the  semi-major  axis 
and  right  ascension  are  observed  in  Figure  23.  Generally,  a  least  squares  method  is  applied  to 
correct  the  initial  conditions  to  remove  these  biases.  Here  we  apply  the  short  period  corrections 
of  Eq.  (1 1 1)  at  the  initial  time  to  remove  the  biases,  as  shown  in  Figure  24. 
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Figure  23.  Comparisons  Between  Numerical  and  Analytical 
Simulations  Without  Corrections  for  the  Initial  Conditions  of 
the  Analytical  Model 
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Figure  24.  Comparisons  Between  Numerical  and  Analytical 
Simulations  With  Corrections  to  the  Initial  Conditions  of  the 

Analytical  Model 


4.2.2.2  Relative  Errors  for  Formation  of  Identical  Satellites 

The  deputy  satellite  has  the  same  area  and  mass  as  the  chief  satellite.  The  mean  orbit  of  the 
deputy  is  determined  from  e0  =  e0  +  Ae0.  The  initial  phase  angle  a0  is  selected  to  be  zero  and 
the  relative  orbit  size  p  is  1  km.  Using  the  corrected  initial  conditions,  Eq.  (107)  is  integrated 
for  the  chief  and  deputy  satellite,  respectively.  Then  nonlinear  relative  position  and  velocity  are 
obtained  by  the  unit  sphere  approach  [18],  which  are  the  reference  values  for  evaluating  the 
accuracy  of  the  extend  GA  STM. 

Figures  25  and  26,  respectively,  show  the  effects  of  differential  SRP  due  to  orbital  element 
differences  for  a  formation.  The  relative  position  error  due  to  the  use  of  a  two-body  STM  is  less 
than  0.06  m  in  10  days,  and  it  is  reduced  to  0.03  m  by  modeling  differential  SRP  due  to  orbital 
element  differences.  Notice  part  of  the  relative  errors  is  from  nonlinearity.  The  theory  does  not 
include  the  effects  due  to  differences  in  solar  areas  of  the  satellites,  which  is  a  more  significant 
effect. 
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Figure  25.  Relative  Motion  Errors:  Two-Body  STM  Vs 
Nonlinear  Simulation 


Figure  26.  Relative  Motion  Errors:  STM  Including  Differential 
SRP  Due  to  Orbital  Element  Difference  Vs  Nonlinear 
Simulation 
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4.2.2.3  Relative  Position  Errors  From  Differential  Area 


Now  we  assume  the  cross  sectional  area  of  the  deputy  satellite  is  5.07  m2  or  about  1%  more  than 
that  of  the  chief  satellite,  as  shown  in  Figures  27-28. 


Figure  27.  Relative  Motion  Errors:  STM  Including  Differential 
SRP  Due  to  Orbital  Element  Difference  Vs  Nonlinear 
Simulation 
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Figure  28.  Relative  Motion  Errors:  STM  Including  Differential 
SRP  Due  to  Orbital  Element  And  Differential  Area  Vs 
Nonlinear  Simulation 
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Figures  27  and  28  show  relative  motion  errors  between  the  solutions  of  the  extended  GA  STM 
and  nonlinear  simulations.  The  effect  of  differential  area  is  modeled  into  the  solutions  of  the 
extended  GA  STM  in  Figure  28  while  the  differential  SRP  due  to  only  orbital  element  difference 
is  included  in  the  solutions  of  the  extended  GA  STM.  Comparing  Figure  27  with  Figure  28,  one 
can  see  the  relative  motion  errors  due  to  the  effects  of  differential  area  are  much  larger  than  those 
due  to  orbital  element  differences.  The  reason  is  described  as  follows:  taking  the  normal 
direction  of  the  acceleration  due  to  SRP  as  example  and  differentiating  Eqs.  (105-106),  we  have 


Pk-r 

m 


8Uh  = - Uh8A - A8Uh 


Pk 
m 


Pk  (  -  8  A 

U ^ - h  C^8i  +  C2dkl 

A 


Am 


V 


(168) 


where 


Cj  =  sin  Q  cos  i  cos  Ae  -  cos  Q  cos  i  sin  Ae  cos  s  -  sin  i  sin  A,  sin  s 
C2  =  cos  Q  sin  i  cos  A,  +  sin  Q  sin  i  sin  Ae  cos  s 

From  Eq.  (112),  we  have 


8i  =  Si0  +  8i 


(169) 


The  long  and  secular  rates  are  assumed  constants.  The  estimated  values  are  approximately  10'16/s 
for  the  inclination  and  right  ascension  right  rates  at  the  initial  time.  Substituting  Eq.  (112)  into 
Eq.  (169),  then  Eq.  (168)  after  ignoring  the  rate  terms,  we  have 


8Uh  = 


Pk 


f 


Am 


\ 


V  h  A  1  ao  J 


U, 


(170) 


Equation  (170)  clearly  demonstrates  even  one  percent  of  the  area  variation  will  dominate  the 
differential  SRP. 

Similar  to  the  drag  problems,  modeling  differential  area  into  the  STM  is  very  important  as  this 
can  be  seen  in  the  same  scale  of  Figures  27-28.  The  root  mean  square  (RMS)  is  0.686  m  with 
modeling  differential  area  while  it  is  4.18  m  without  modeling  differential  area  into  the  extended 
GA  STM  for  one  percent  difference  between  the  cross  sectional  areas  of  the  chief  and  deputy 
satellite.  Table  2  shows  the  effect  of  modeling  differential  area  into  the  STM  to  deal  with  the  SRP 
perturbations  due  to  different  differential  area  among  the  satellites  and  presents  the  RMS 
variations  as  the  percentage  of  the  area  difference  increases  for  a  10  day  simulation.  There  is 
about  six  times  reduction  for  relative  motion  errors  for  the  example. 


Approved  for  public  release;  distribution  is  unlimited. 

73 


Table  2  RMS  Variation  With  Area  Difference 


Deputy  Area  (m2) 

Area  Diff  (%) 

RMS  without  Diff 

RMS  with  Diff  Area 

Area  (m) 

(m) 

5.07 

1 

4.18 

0.69 

5.12 

2 

8.36 

1.37 

5.17 

3 

12.54 

2.06 

5.22 

4 

16.71 

2.74 

5.27 

5 

20.90 

3.43 

4.2.3  Perturbations  due  to  Drag 

The  Chief  satellite  orbit  is  a=6700km,  e=0.004,  i=48deg,  Q=0,  eo=10  deg,  and  M=\20  deg.  The 
simulation  time  is  3  days.  The  satellite  is  100  kg  and  the  cross  sectional  area  is  1  m2.  The  relative 
orbit  size  is  1  km  and  initial  phase  angle  is  zero.  The  deputy  satellite  orbit  is  determined  in  the 
usual  manner  for  the  PCO.  All  the  initial  conditions  are  obtained  by  analytical  methods. 

4. 2. 3. 1  Drag  Only  Problem 

We  use  the  two-body  GA  STM  to  get  the  solutions  of  relative  motion,  and  integrate  nonlinear 
equations  including  drag  perturbation  only  for  the  chief  and  deputy  orbits  and  get  the  real 
relative  solutions.  The  errors  are  the  differences  between  the  former  and  latter,  as  shown  in  Fig. 
29.  Then  we  analytically  model  the  drag  into  the  extended  GA  STM,  the  corresponding  errors  are 
shown  in  Fig.  30. 


Figure  29.  Relative  Position  Errors  From  Two-Body  Modeling 
For  Drag  Only  Problem 
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Figure  30.  Relative  Position  Errors  From  Modeling  Drag  Into 
Extended  GA  STM  for  Drag  Only  Problem 


There  is  more  than  one  order  of  magnitude  improvement  by  including  drag  in  the  extended  GA 
STM,  as  shown  in  Figs.  29-30.  Since  the  assumption  on  the  constant  rates  perturbed  by  drag  is 
inaccurate  [44],  the  accuracy  is  not  high. 

The  cross  sectional  areas  are  the  same  for  the  chief  and  deputy.  Now  we  assume  there  is  one 
percent  difference  in  the  area,  i.e.  the  area  of  the  deputy  is  1. Ole-6  m2  and  the  chief  area  is  still 
1.00e-6  m2.  The  simulation  time  is  one  day,  as  shown  in  Figs.  31-32. 
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Days 


Figure  31.  Relative  Position  Errors  Without  Modeling  Diff. 
Area  Into  Extended  GA  STM 
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□  ays 


□ays 


Figure  32.  Relative  Position  Errors  With  Modeling  Diff.  Area 
Into  Extended  GA  STM 


Figures  31-32  indicate  the  differential  area  has  significant  effects  on  the  accuracy,  especially  in 
the  in-track  direction.  Modeling  the  differential  area  into  the  extended  GA  STM  remarkably 
alleviates  the  effect.  As  the  differential  area  percent  increases,  it  is  necessary  to  modeling  the 
differential  area  into  the  extended  GA  STM,  as  shown  in  Fig.  33  for  a  one  day  simulation. 


□  iff  Area  Percent 


Figure  33.  Effect  of  Modeling  Differential  Area  in  the  Extended 

GA  STM 


Figure  33  demonstrates  it  is  important  to  introduce  the  area-mass-ratio  variable  Ct  into  the 
extended  GA  STM,  especially  for  the  perturbations  sensitive  to  the  cross  sectional  area. 
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4. 2. 3. 2  Drag  and  Gravity  Combined  Perturbations 

Comparing  Figs.  29-30  with  Figs.  34-35,  one  can  see  the  effect  of  the  gravity  perturbation  is 
about  four  times  more  than  that  of  the  drag  perturbation  for  this  example.  The  accuracy  modeling 
J2  is  much  higher  than  that  modeling  the  drag  perturbations  since  the  mean  rates  caused  by  the 
gravity  perturbations  are  exactly  constants. 

The  accuracy  of  the  extended  GA  STM  is  shown  in  Fig.  36  for  the  drag  and  J2  combined 
perturbations. 


□ays  Days 


Figure  34.  Relative  Position  Errors  From  Two-Body  Modeling 
For  J2  Only  Problem 
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Figure  35.  Relative  Position  Errors  From  GA  STM  For  J2  Only 

Problem 
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□ays 


□ays 


Figure  36.  Relative  Position  Errors  From  Modeling  Combined 
Perturbations  Without  Height  Corrections 


The  errors  for  the  combined  perturbations  are  larger  than  those  obtained  by  a  superposition  of  the 
individual  effects,  due  to  a  coupling  effect. 


□ays 


□ays 


Figure  37.  Relative  Position  Errors  From  Modeling  Combined 
Perturbations  With  Height  Corrections 


Fig.  37  illustrates  the  coupling  effects  are  greatly  reduced  by  applying  the  height  correction  given 
by  Eq.  (132).  The  cross  sectional  areas  are  the  same  for  both  the  chief  and  deputy  for  the 
combined  perturbations.  Now  we  want  to  see  the  effect  of  the  differential  area  on  the  drag  and 
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gravity  combined  perturbations  for  one  day  simulation.  Fig.  37  shows  the  RMS  position  errors 
vary  with  the  differential  area  percent  changes. 


□  iff  Area  Percent 


Figure  38.  Effect  of  Modeling  Differential  Area  Into  Extended 
GA  STM  For  Combined  Perturbations 


□ays  Days 


Figure  39.  Relative  Position  Errors  From  Modeling  Drag  and 
J2  By  Semi-Analytical  Method 


Fig.  38  indicates  there  are  little  changes  if  the  combined  perturbations  are  considered  for  the 
effect  of  the  cross  sectional  area,  possibly  due  to  gravity  perturbations  being  insensitive  to 
variation  of  the  area.  Figure  39  shows  the  results  of  using  the  semi-analytic  STM.  Comparing 
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with  Fig.  38,  one  can  see  the  relative  position  errors  are  effectively  reduced  by  using  the  semi- 
analytic  STM. 

4. 2. 3. 3  Effects  of  Different  Density  Models 

We  only  consider  perturbations  from  drag  only.  All  the  simulations  are  nonlinear.  We  integrate 
the  equations  of  motion  for  the  drag  only  problem  for  chief  and  deputy,  respectively,  and  then  get 
relative  distances  that  are  compared  with  the  GMAT  solutions  with  the  same  initial  conditions. 
The  density  model  is  the  Jacchia-Roberts  model  in  GMAT  simulations.  We  use  two  kinds  of  the 
density  models,  the  first  one  is  the  1976  US  standard  atmosphere  and  the  second  is  Harris- 
Priester  density  model.  [45] 

The  epoch  time  for  GMAT  simulations  is  April  15,  00:00:00,  2015.  Fig.  40  shows  the  relative 
distance  error  from  nonlinear  two-body  simulation  and  GMAT  simulation  with  Jacchia-Roberts 
model.  Fig.  41  illustrates  the  relative  distance  error  from  nonlinear  simulation  with  1976  US 
standard  density  model  and  GMAT  simulation  with  the  Jacchia-Roberts  model.  Fig.  42  shows  the 
relative  distance  error  from  nonlinear  simulation  with  the  Harris-Priester  density  model  and 
GMAT  simulation  with  Jacchia-Roberts  model. 

The  standard  density  model  only  considers  a  spherical  atmosphere  while  the  Harris-Priester 
density  model  introduces  the  effects  of  oblateness  and  diurnal  bulge.  From  Figs.  40-42  we  can 
see  the  nonlinear  modeling  errors  obviously  decrease  as  we  use  more  accurate  density  models. 
Assuming  we  exactly  model  nonlinear  systems  by  the  extended  GA  STM  with  the  standard 
density  model  or  Eq.  (116),  there  is  about  an  0.8  km  errors  at  the  end  of  three  day  simulation  if 
we  compare  with  Jacchia-Roberts  model,  as  shown  in  Fig.  41. 


Figure  40.  Relative  Position  Errors  By  Using  Two-Body  and 
Jacchia-Roberts  Model 
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Figure  41.  Relative  Position  Errors  By  Using  Standard  Density 
and  Jacchia-Roberts  Model 
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Figure  42.  Relative  Position  Errors  by  Using  Harris-Priester 
and  Jacchia-Roberts  Model 
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4.3 


Navigation  and  Thruster  Inaccuracies 


Figure  43.  Time  history  of  -  1  (lines)  and  m'  -  1 

(crosses)  color  coded  by  i  (1  =  blue,  2  =  green,  3  =  red,  4  =  teal, 
5  =  purple,  6  =  gold)  for  all  j.  Values  plotted  once  every  orbit 

for  clarity. 


As  a  practical  example,  suppose  the  true  dynamics  include  zonal  terms  of  the  Earth  gravity  up  to 
sixth  order  whereas  the  approximate  dynamics  include  those  only  up  to  second  order.  Figure  43 
is  the  time  history  of  O'T/O'7  -  1  compared  with  m'  -  1  over  40  orbital  periods  for  an  object 
whose  orbital  elements  are 

(a,  e,  i,  Q,  a>,  M) 

=  (7554.9  km,  0.050078, 0.83775  rad,  0.34859  rad,  0.1 7984  rad,  2.089  rad) 

at  epoch.  We  see  that  the  two  plots  match  well  for  the  entire  analysis  timespan  of  3.0256  days. 

We  further  find  in  Fig.  43  that  the  difference  between  values  in  the  two  STMs  are  on  the  order  of 
0.1  %,  suggesting  that  a  linear  function  relating  STM  bias  and  maximum  a  posteriori  state 
estimates  will  sufficiently  capture  sensitivities  between  them  for  many  relative  navigation 
scenarios  of  interest. 
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We  now  demonstrate  how  the  tools  developed  in  this  research  enable  rapid  design  of  relative 
navigation  algorithms.  Consider  a  formation  of  two  spacecraft  where  the  initial  orbital  elements 
for  the  chief  are  as  given  and  the  differences  in  orbital  elements  between  the  chief  and  deputy  at 
epoch  are 

(A  a,  Ae,  A  i,  AO,  Aoj,  AM) 

=  (0.00055381  km,  5.7624  x  10"5, 1.0515  x  10“5  rad, 

-  9.2463  x  1(T8  rad,  -4.939  x  1(T6  rad,  5.0279  x  1(T6  rad). 

For  this  work,  we  assume  the  absolute  state  of  the  chief  is  perfectly  known  and  that  the  true  and 
reference  relative  trajectories  are  set  equal.  This  assumption  shall  be  relaxed  in  future  work. 

The  observations  are  instantaneous  range  measurements  made  every  0.081001  seconds, 
corresponding  to  approximately  100  measurements  over  the  course  of  a  single  orbital  period  for 
the  chief.  Measurement  noise  is  assumed  to  be  Gaussian  white  noise  with  a  constant  standard 
deviation  of  1  meter.  The  measurements  are  simulated  over  3.0256  days  or  approximately  40 
orbital  periods  of  the  chief,  resulting  in  4001  range  measurements  total. 

A  batch  processor  is  run  for  the  true  case  whose  dynamical  model  includes  a  40  x  40  gravity 
field.  State  estimate  and  covariance  biases  are  then  computed  for  three  approximate  models 
which  includes  a  20  x  20,  10  x  10,  or  2  x  0  gravity  field.  Any  deviation  from  the  estimated 
trajectory  based  on  the  true  dynamics  would  indicate  extraneous  fuel  expenditure  attributed  to 
dynamical  model  error.  The  total  deviation  AX  is  the  sum  of  the  state  estimate  bias  and  the 
deviation  in  the  reference  state 


AX(r)  =  Ax(r)  +  Ax(0  (171) 

This  vector  may  be  converted  into  an  approximate  Av  cost  as  a  consequence  of  Gauss'  planetary 
equations.  Suppose  AX(t)  =  (8a,  8e,  Si,  8(1,  Son,  8M )  are  expressed  in  the  mean  orbital 
elements.  An  osculating-to-mean  transformation  including  first-order  ,h  effects  such  as  Brouwer 
theory  [3]  will  suffice  for  this  computation  regardless  of  dynamical  model  used.  Then,  assuming 
that  the  sensitivity  between  the  mean  and  osculating  elements  are  sufficiently  close  to  unity, 


A Vf,  =  (h/r)  V Si2  +  6d2  sin2  i 
A vrp  --  ~(naj 4) 

A  Vra  =  ~(naj  4) 


(1  +e)2/rj 

(\-eflrj 


(■ 8a>  +  cos  i)  +  SMj 
(6u>  +  <5Q  cos  i)  +  8M\ 


Avtp  =  nar}IA(5aja  +  5ej{  1  +  <?)] 
Av/o  =  nar](A  [Saja  -  Se/(  \  -  e)] 


(172) 


where  subscripts  h,  r,  and  t  indicate  that  the  thrust  is  made  along  the  angular  momentum 
direction,  radial  direction,  or  a  direction  perpendicular  to  the  two,  respectively,  and  subsubscripts 
p  and  a  indicate  that  the  maneuver  is  executed  at  periapsis  or  apoapsis,  respectively. 
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The  state  estimate  uncertainty,  however,  may  be  so  high  such  that  AX  cannot  be  distinguished 
from  noise.  Here,  noise  is  attributed  solely  to  measurement  error,  but  in  subsequent  studies,  we 
hope  to  add  the  effects  of  process  noise.  We  would,  nonetheless,  like  to  compute  another  metric 
characterizing  the  magnitude  of  AX  with  respect  to  true  state  deviation  estimate  x;  i.e.,  if  AX  is 
much  smaller  in  magnitude  compared  to  representative  values  of  x,  then  it  may  safely  be 
ignored.  We  define  a  non-dimensional  distance  d,  similar  to  the  Mahalanobis  distance  [49], 
where  AX  is  normalized  by  the  diagonalized  true  state  estimate  covariance  F 


d(t )  = 


AX(t)TP(t)-1AA'(t) 


(173) 


Thus,  a  d  metric  value  of  1  would  indicate  that  the  magnitude  of  A X(t)  in  each  coordinate 
direction  is  equal  to  their  respective  true  l-cr  uncertainty  level.  We  choose  to  ignore  the 
correlations  in  P  as  we  are  not  interested  in  quantifying  the  direction  of  AX(t)  with  respect  to  the 
state  uncertainty  for  this  particular  metric. 

We  find  for  this  problem  that,  in  under  10  orbits,  A Hi  becomes  small  such  that 


A xvPrj^P^HjRt'Ay,]) 


(174) 


based  on  Eqs.  (160)  and  (164).  17  is  nearly  parallel  to  N  =  £  HjR^yi  as  the  STM  term 
dominates  (4).  Therefore,  after  running  the  true  batch  processor,  one  only  needs  to  compute  Ayi, 
which  is  often  available  as  a  closed- form  function  of  the  solution  flow,  and  subsequently  17 
across  approximate  dynamical  models  to  evaluate 

AX(t )  «  <E>(t;  Xq)Pi 7  +  Ad>(t;  x0)x0  (175) 

accurate  to  order  of  magnitude.  Table  3  corroborates  this  reasoning  based  on  numerical 
simulations.  That  is,  the  truncated  expression  Eq.  (175)  for  A X(t)  results  in  the  d  distance  metric 
computed  accurately  to  at  least  order  of  magnitude. 

Table  3  compares  for  each  approximate  dynamical  model,  the  non-dimensional  distance  metric  d 
using  the  truncated  form  of  AX(t)  in  Eq.  (175)  (t/irunc)  and  subsequently  substituting  into  Eq. 
(173),  and  by  running  individual  batch  processors  to  compute  AX(t)  (<7ruii).  A  is  the  relative 
error  of  the  two  methods  in  terms  of  percentage. 
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Table  3  Metric  Comparison  For  Various  Gravity  Models 


Gravity 

^Trunc 

*/fu11 

A  [%1 

20  *20 

5.0496  x  10  2 

4.6012  x  lO'2 

9.7454 

10  x  10 

2.8434  x  10'1 

2.8231  x  lO’1 

0.7177 

2*0 

2.5473  x  10° 

5.7479  x  10° 

-55.682 

For  both  the  20  x  20  and  10*10  models,  the  total  deviation  of  the  state  estimate  due  to 
inaccurate  dynamical  modeling  is  an  order  of  magnitude  smaller  than  the  magnitude  of  the  state 
estimate  noise.  If  we  are  to  implement  these  in  our  navigation  algorithms,  then,  it  is  unlikely  that 
fuel  will  be  spent  to  correct  effects  due  to  unmodeled  terms  in  the  geopotential.  The  same  cannot 
be  said,  however,  for  the  2  *  0  model,  as  the  magnitude  of  the  bias  in  the  state  estimate  is  well 
over  the  3-cr  level  of  the  diagonalized  state  covariance,  as  indicated  by  c/ruii  =  5.7479.  For  the 
given  measurement  parameters,  then,  a  2  *  0  model  may  be  deemed  insufficient  as  it  leads  to 
wasted  fuel.  On  the  other  hand,  the  dimensional  value  of  AX(t )  throughout  the  analysis 
timeframe  is,  at  most,  meter-level  in  position  and  mm/s-level  in  velocity.  Qualitatively,  we 
expect  that  correcting  such  a  small  state  deviation  is  within  the  bounds  of  thruster  error.  Future 
work  is  to  directly  relate  the  Av  value  based  on  Eq.  (172)  and  control  system  errors. 

Finally,  we  note  that,  if  the  sensor  were  precise  to  0. 1  meters  1-c,  then  without  running  any 
additional  simulations,  linearity  suggests  that  d  will  accordingly  increase  by  approximately  an 
order  of  magnitude  compared  to  those  in  Table  3,  bringing  the  bias  of  the  10  *  10  model  close  to 
3-cr  level  of  the  estimate  noise. 

We  have  thus  shown  how,  with  the  tools  developed  in  this  work,  one  is  able  to  bypass  running 
multiple  batch  processors  per  each  dynamical  model  we'd  like  to  test  in  navigation  algorithm 
design.  Further  computational  gains  may  be  reaped  by  applying  analytical  STMs  or  by 
computing  the  STM  bias  based  on  a  vector  scaling  function  between  said  models.  In  order  to 
apply  this  methodology  to  a  wider  range  of  measurement  types  and  mission  scenarios,  we  are 
working  to  develop  a  way  to  automatically  determine  ignorable  terms  in  the  linear  sensitivity 
function. 
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5  CONCLUSIONS 


The  conclusions  of  this  program  sponsored  by  the  Air  Force  Research  Laboratory  are: 

•  A  systems  approach  for  the  selection  of  the  dynamic  model  that  is  based  on  the  accuracy 
of  the  relative  navigation  system  and  the  accuracy  of  the  thrusters  has  been  developed. 
The  methodology  simplifies  the  workflow  of  selecting  the  dynamic  model  and  designing 
navigation  systems  so  that  the  trade  space  between  navigation  system  parameters  and 
dynamical  model  fidelity  could  be  quickly  surveyed  in  lieu  of  performing  massive 
numerical  simulations  for  numerous  scenarios  and  system  parameter  variations. 

•  Including  the  short  period  terms  of  higher  order  geopotential  effects,  especially  J\  ,  in 
calculating  the  relative  state  initial  conditions  can  reduce  the  secular  error  drift. 

•  Including  the  secular  and  long  period  effects  of  the  lunar  perturbations  can  improve  the 
accuracy  of  the  relative  state  prediction  for  high  altitude  formations. 

•  The  correction  of  the  deputy  semi-major  axis  for  negating  in-track  secular  drift  has  been 
expanded  to  include  higher  order  geopotential  terms  and  this  expansion  will  further 
reduce  the  drift. 

•  The  GA-STM  uses  a  set  of  nonsingular  variables  for  non-equatorial  orbits  and  equinoctial 
elements  for  near  equatorial  elements.  Other  sets  of  nonsingular  variables,  such  as  Hoots 
variables,  were  evaluated  and  compared  to  those  in  the  GA-STM.  In  some  cases  the 
performance  was  better  and  others  it  was  not.  No  general  methodology  for  evaluating 
different  sets  of  variables  was  found. 

The  recommendations  from  this  project  are: 

•  The  system  methodology  for  selecting  the  dynamic  model  and  relative  navigations 
system  based  on  relative  navigation  requirements  and  thruster  accuracy  should  be  used  to 
reduce  the  analysis  required  in  making  these  selections  in  formation  design. 

•  To  help  reduce  secular  drift  the  higher  order  geopotential  short  period  terms  should  be 
included  in  the  calculation  of  the  relative  state  initial  conditions. 

•  The  extensive  analysis  performed  in  this  project  has  revealed  the  need  for  Air  Force 
Research  Laboratory  to  have  a  detailed  formation  system  simulation  that  includes  all 
relevant  perturbing  forces,  control  forces  and  modules  for  thruster  control,  and  modules 
for  different  types  of  relative  navigation. 
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APPENDIX  A.  FORMULAE  FOR  ZONAL  HARMONICS 

A.l  Expansion  Formulae 


The  following  formulae  are  helpful  in  computing  closed-form  expressions  for  the  integrals  given 
in  Eq.  (6).  These  formulae  include  the  definition  of  Legendre  polynomials  and  formulae  for 
converting  powers  of  trigonometric  functions  into  sums  of  their  arguments. 
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A.2  Intermediate  Terms  for  Short-Period  Corrections 


The  following  intermediate  terms  are  used  in  the  analytic  formulae  for  short-period  contributions 
due  to  an  arbitrary  zonal  harmonic  with  degree  greater  than  two. 
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=  2^  (• lcos(/> + ecos(2/> + e2,)  i++^+  3)  *2-» 


LVl-J 


LfJ  n—l  ,  v  , 

<-  /  afc-l,n(e)  W  —  K 

Ma’e)  /  ,Pj,n(l)  /  „  Ok  h 

j=0 


even[odd]  /c=l 


2fc  A; 


(A.  14) 


/  fc  \ 

y:  7pJ,n  (  k—{n—2j—2p)  )  (n  -  2J  -  2p)(~  sin[cos]  (n  -  2j  -  2p)g) 

P= 0  V  2  / 


Lfj  r  j  l^j-j 

=  £(a,e) £&,„(»)  j  (4cos(/)  4-  ecos(2/)  +  3e)  £  7Pj,n  cos  [sin]  (c/  +  eg) 


j=0 


p= 0 


Qfc-I,n(e)  n  -  k  rf  f  dW3  df  _  1  dWA  \  1 
^  2k  k  y/Jui  \  df  dl  g  dg  )  J 


fc=l 


(A.  15) 
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APPENDIX  B.  O  MATRIX  FOR  THIRD  BODY  PERTURBATIONS 


Derivatives  for  a 


Derivatives  for  6 1 


da 

8a0 


da  _  da 
di0  do0 


ao0 


(B.l) 


d0_ 

da0 

de_ 

de0 

de__ 

h= 

de 

Ho 

de 

Ho 

de 

dQ,, 


dX_ 

da0v  '0/  dqx  da, 


('-0- 


dG  dqx  dG  dq2 


o 


H  H 


dG 

de 


dG 

v  H  j 


dG 

de 


dX  /  \  dG  dq 

r'-'o- 


dG  dq 


di, 


0 

di 

Ho 

di 

dq 


H  H  H  H 

dG  dqx  dG  dq 

H  Ho 

dG  dq. 


dG 


de 


H  Ho 

dG  dq2 


20 


H  Ho  H  H 


20 


J 

rdG_" 

Kde  j 

rdG^ 

de 


(B.2) 
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Derivatives  for  i 


di  di 


da0  da0 

3i  =  0 


d0n 


di  .  di  /  \ 


0 

di 


A 

di 


Al0  Al0 

di  di 

^#20  ^20 

di  di 

d£2  d£2 

w  o  u  o 


('-o 

('-o 


Derivatives  for  qi 


(A 

d  — 

A  =  KK  —+k  A£Z 

da0  1  3da0  2  da0 
dq 


den 


L-0 


'e' 


dq.  d(b 

—  =  k]k3—+ k2 


\eJ 


fe' 


H  =  K  K  dCQ 

A10  Al0  A 

dq,  dd) 

-^  =  KtK- —  +  K7  ^ 
dq2()  dqw  dq 


\eJ 


10 

fe' 


vey 


+  K,  cos  Aft) 


K}  sin  Aft) 


20 

<bL=KK<M+KA!U 

dQ0  3  dQ0  2  dQ0 


(B.3) 


(B-4) 
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Derivatives  for  qi 


Derivatives  for  Q 


H  =  K  K  3®. 

3a„  2  3  3a„ 


vgy 

1  3a„ 


3g 


^  =  0 


3e„ 


3a,,  3® 

' dL=KlK3dL~K ' 


1  3/„ 


a^2 


3ft) 


v^y 


Ho 

2  3  Ho  *  Ho 

(  -\ 

dq2 

3 

e 

=  K^K,  d°J  K,  . 

UJ 

H0 

^=K  K _ 

3Q„  2  3  3Q 


5?20 

3g) 


^20 

( 

3 


+  AT,  sinAto 


+  Ay  cosAtw 


K, 


\eJ 

3Q„ 


(B.5) 


3Q  3Q 

3a0  3a0 

^  =  0 
30n 


(»-<.) 


3Q  3Q  /  \ 

3r*(,-° 

3Q  3Q 


a^io  3?io 

3Q  3Q 


^20  ^20 

3Q  „  3Q 


3Qn 


1  +  - 


('-o 

('-'.) 


3Q 


(B.6) 
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where 


G-l-k[t0)~  M[tQ)  +  (d{tQ)  +  (M  +  (i))(t-tG) 
K]  =(-glosinA£0  +  ^2OcosA(w) 

K2  =  (g10cosAfy  +  g20  sinAto) 


(B-7) 
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APPENDIX  C.  O  AND  D  MATRIX  FOR  SRP  PERTURBATIONS 
C.l  O  Matrix 


Derivatives  of  a 


Derivatives  of  9 


a  =  aQ 

da 

Sa0 

da  da 

da 

da  da 

da 

n 

d00  di0 

dch. 

5q20  50 Q 

'dCB~ 

59 

(  52 

—  | 

5G  5qx 

5G  5q /, 

5G 

da0 

[5a0 

5ql  5a0 

dq2  daJ/ 

59 

59 

5G 

/  5G 

'  59 

59 

(  52 

—  | 

5G  5qx 

5G5q. ^  jdG 

di0 

\di0 

5qx  5i() 

dq2  5i0 )  j  59 

59 

f  52 

5G 

5G  5qx  5G 

dq2] 

jdG 

d%< 

dqw 

5qx  5qX0  5q2 

dqj/ 

/  59 

59 

(  52 

5G 

5G  5qx  5G 

dqA 

/dG 

d(h0 

^1^20 

dqw 

5qx  dq2Q  dq2 

dq2Ji 

/  59 

59 

f  52 

_  —  | _ 

5G  5qx 

5G  dq2  \ 

/dG 

5qx  50  ( 

,  dq2  5QJ/ 

59 

59 

f  52 

—  | 

5G  5qx  5G  dq2  \ 

/dG 

dCB 

6qx  5Cb  dq2  5Cb)i 

/  59 

(C.l) 


(C.2) 
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Derivatives  of  i 


Derivatives  of  qq 


1  ~  h  +  hA1 
di  3PkCB  -  A 
da0  ~  4 ml  ' 


di  ,  .  3PkCBq .  dUh 


— -  =  — - — PkC  U  At 
dqw  2na0  u  “ 

f=0 

dcho 

di  =  3PkCBq{()  ?)Uh 
5Q0  2  naQ  dQ0 

di  3 Pk  - 
dCB~  2na0  A  1 


^  =  ^io  +  ^ioAf 

H  _  3  , 


PkCJJAt 


— —  PkC  ^-L  At 
2nan  dL 


- PkC  — — -  At 

2nan  B  dQn 


(C.3) 


(C.4) 
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Derivatives  of  q2 


Derivatives  of  Q 


42=  420+^20^ 

PkCBU  At 

dan  Anal 


dq  3  ,  dU 

-^T  =  ~ - PkC  At 

dL  2na„  dL 


3Q0  2 naQ  B  3Q0 


3<?2  _  3 Pk 
dCB  2na0 


£2  —  £20  +  Q(JA  t 

= - y - PkCUqAt 

da0  Anal  sin  i()  "  '  ' 


3£2  3  PkCBq2Q  -  .  at/. 

—  = - gm-  C/.coti iL 

oz0  2«a0sin/0^  o/(| 


= - PkCUAt 

dq20  2naQ  sm  i()  '  " 

dQ.  3  3t7A 

— —  =  1 H - PkCRq7,  — — - 

3£2„  2na„sinz'  0  3£2„ 


3C_  2«flnsin/n 

n  U  U 


Uk4  2At 


(C.5) 


(C.6) 
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where 


dX  3 n  3PkC  /  -  -  -  \ 

aX ^T{  ’g'° +  + tW°  ^ 


ex  3PkCn  ( dU  8U 


dL  2  na^  l  di 


f*  10  +  ?20  cot/n  I  Ar 


0  V  ^‘0 


y20  wl,o 

Ul0 


PkCJJ  At 


dqw  2na 


dq7,  2na, 


aQ„  2nan  H  aQ„  10  20  aQ„  20  °, 


■«C„((7,  +  C/,coti0)A( 

,  ^  f  a£7  dut  du,  } 

-PkCj - -q. nH - Lo,j - -qincotL  At 

B\  ^10  ^20  ^20  0 


dX  3  Pk 

dCB  2  naQ 


(UrCh0+U,Ch0+Uk(l20COti0)At 


(C.7) 


Derivatives  of  Cb 


ac,  =  ac,  =  ac,  =  ac,  =  ac,=0 

a«0  5/0  57l0  ^20 


(C.8) 
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C.2  D  Matrix 


Derivatives  of  da 


^  PkC 

2  ^(-sinl  -  qx  sin 2 1  +  q2  cos 2 A^Ur  +  (cosl  +  q2  sin 2 1  +  qx  cos 2/1  J 


^=- 


0£>a 

5a 

0£a 


6PkCr 


n2  a 


cosl  +  ~q±  cos2A  +  ^q2  sin21  j  Ur  + 


sin  A  -  ^  g2  cos  21  +  ^  qx  sin  21 


[/ 


01 

- =  a,  — 

0£  2  0# 


0£a 

di 

88  a 
dqx 

88  a 
dq2 

dSa 

~dQ. 

dSa 

dcT 


2PkC„ 


n 


1  1 

cos  A  +  —  qx  cos  21  + — q2  sin  21 

v  2  2 

01 


N  at/ 


0/ 


-  + 


sinl-^2  cos  21  ~i~  ~qx  sin2lj  -^-x 


P£CR  /  -  -  X 

- -±(Ur  cos 21  +  [/  sin 21) 

n  v 


-ha. 


PkC„  i  -  -  \ 

B-\Ur  sin21  - U t  cos21j  +  < 


n 


8qx 

01 

dq2 


2PkC, 


cosl  ~h  ~q^  cos21  +  ~q2  sin21 


8U 


2Pk 


cosl  ~h  ~q^  cos21  +  ^2  sin21 


)  0Q 

r 


-  + 


sin  1  -  ^  q2  cos  2 1  +  ^  qx  sin  2 1 


\  CST  T 


dU. 


dQ 


U  + 

r 


1 


1 


sin /1-— q2  cos  2/1  +—^  sin2/l 


U. 


(C.9) 
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Derivatives  of  86 


dse 

da 

886 


U  ,  O  0  ,  o  •  AdSA  1  O  0 

=  1  +  2 q  cos  A  +  2q  sin  A  - +2sin/l - -  -  2cos/l 

v  '  da  da 


da  da  da 

(.  .  ,  .  .\88A  .  .  ndSq.  ~d8q 

=  11  +  2 q.  cos  A  +  2 q,  sin  A I - +2smA - -  -  IcosA - - 

86  1  ’  86  86  86 

/  .  .\3/l  _  ,  8 A  „  .  .  ,8A  „ 

+2  -(7,  sin  +  +  (7,  cos  /l  —  cm,+2cos/1 — <?<7,  +  2sin/l — do, 

U  ,  o  o  ,  o  •  ,  dSq2 

=  1  +  2 q.  cos  /l  +  2 <7,  sin  /l  - +2sin/l - 1  -  2cos/l - - 

v  'da  di  di 

(.  .  ,  „  .  x  5t)2  .  .98q 

=  11  +  2^  cosA  +  2q  sm/l  - +2sm A - 

v  '  dqx  dq 


886 

di 

886 


dqx 

886 

dq2 

986 

80. 

886 

8CD 


1  -  2cos/l  — — +  28A  cos  A 


8q. 


=  ( 1  +  2 q.  cos /l +  2(7,  sin/l)^^+2sin/l-— ^ -AcosA^—^1  +28AsinA 
V  yi  >  9a  9a  9a 


8q 

88A 


dq2 

98 q. 


dq2 

80 


=  ( 1  +  2 q.  cos  A  +  2 q,  sin  A)  -  +2sin/t  _  2cosA 

1  yi  ’  80  80 

.  ,  \  88A  .  98 q.  88q . 

=  1  +  2(7,  cos  A  +  2(7,  sm  A  - +2sin/l - -  -  2cos/l - : 

1  1  y2  'SCD  8Cd  dCD 


(C.10) 
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Derivatives  of  Si 


h=~ 


PkCU 

_ B  h 


2  n2  a 

dSi__  PkCJU, 
da 


(2 cos  A  +  q{  cos  2 A  +  q2  sin 2/t) 
1  .  _  1 


A 


2  2 

n  a 


s-h-'  2  sin  A  +  ^  qx  sin  2A  cos  2A 


/ 


dSi  _  .  5A 
~h~80 


50 

dSi _ ^ 

di  2  n2  a 


PkCn 


1 


1 


2sin/l  +  —  q  sin2/l  —  q  cos2A 


v 


2  i\  2  2 


at/. 

_ A 

5/ 


a*-  PIC.t/,  .  .  5A 

- = - 7—^  sin  2/1  +  /' — 

4//2a  A  ‘ 


dq. 


5Si  PkC  U,  „ ,  . 

- = - f— ^-cos2/l  +  /. 

4n2a  1 


dq2 

dSi  PkC 

5CI  2 n2 a  K 

dSi 


dql 

5A 

dq2 


1 


1 


A 


5U, 


— 1  2svc\A-\--^qx sin2/i  —  cos2A 

1  .  _  1 


J 


5n 


dCB  2  n  a  ^ 


PkUb-( 2 si nA  +  —  qx  sin 2A  ~—q2 cos 2A 

2  2  ) 


CC.11) 
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Derivatives  of  8qx 


<hx 


dSqx 

da 


PkC  r  - 

— -  -  -  [(2  sin  2A  +  qx  sin  A  +  5 q2  cos  A  - 3 qx  sin  3  A  +  3 q2  cos  3/1)  U r 

(2  cos  2 A  +  3q2  sin  A  -  3 qx  cos  A  +  3q2  sin  3 A  +  3 qx  cos  3 A)  Ut  ] 

PkC 'g 
2  n2a2 


[(-cos 2/1 -g,  cosA  +  5q2  sin/l  +  g,  cos3A  +  q2  sm3A)Ur 


+(sin2/l-3t/2  cos /l -3^  smA-q2  cos 3 A  +  qx  sin 3 X)Ut~^  +  q2  cos i 


den 

da 


dSqx  dA  .  ddkl 

- L  =  ^7ij - HO,  cosz - 

86  n  dO  2  dd 

dSqx  PkCB  r,  .  >  dUr 

- -  = - r-5-  1 -  cos  2  A  -  qx  cos  A  +  5q2  sm  A  +  qx  cos  3  A  +  q2  sin  3  A )  — 51 

di  4  n  a  L  di 

dU 

+(s,m.2A-3q2  cosA-3qx  smA-q2  cosS/l  +  z/j  sin 3 A) — L 

di 

.SO  . 

+q2  cos z - q2  sm loLl 

di 


U  ^ 

dqx 

dA 

- h 

dq  2 


dSqx 

PkCB 

dqx 

4  n2  a 

dSqx 

PkCB 

dq2 

4n2  a 

dSqx 

PkCB 

dQ. 

4zz2a 

+(sin  2  A 

88qx 

Pk  r 

.8X1 

1  dqx 

8X1 

- h 

dq2 


du^ 

8Q 


+(sin2/l-3^2  cos/l-3^  smA-q2  cos,3A  +  qx  sin 3/1) 


dU, 

dQ 


+  q2  cos  i 


.8X1 

dQ 


da 


4  na 


[(-  cos  2 A  -  qx  cos  A  +  5 q2  sin  A  +  qx  cos  3 A  +  q2  sin  3 A )Ur 


+ 


(sin  2  A  -  3q2  cos  A  -  3  qx  sin  A  -  q2  cos  3  A  +  qx  sin  3/1)  Ut~J  +  q2  cos  i 


8X1 

~dC~« 


(C.12) 
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Derivatives  of  Sq2 


<hx 


+ 


ddq2 

da 


PkC  r  - 

-  -  [(2 cos 2A  +  3q2  sin/l-3 ql  cos/l  +  3 q2  sm3A  +  3qx  cos3/l )Ur 

(2 sin 2/1 -5 qx  sin /l -q2  cos/l  +  3^  sin 3/1 -3 q2  cos 3/1  )t/,] 

PkCB 
2  n2a2 


[(sin  2/1  -  3q2  cos  A  -  3 qx  sin  A  -  q2  cos  3 A  +  qx  sin  3 A)  Ur 


+ 


(-  cos  2 A  +  5qx  cos  A  -  q2  sin  A  -  qx  cos  3  A  -  q2  sin  3/1)  Ut  ]  -  qx  cos  z 


oXl 

da 


d5q2  dA  .  ddkl 

- -  =  <793 - <7,  cosz - 

d6  1X  dO  1  dO 


dSq2 

di 


dA 
de 
PkC \ 
4  n2a 


[(sin  2  A  -  3q2  cos  A  -  3 qx  sin  A  -  q2  cos  3 A  +  qx  sin  3  A) 


dUr 

di 


dU 

+  (-  cos  2 A  +  5 qx  cos  A  -  q2  sin  A  -  qx  cos  3  A  -  q2  sin  3 A)  — L 

di 

.dm  . 

-qx  cos  i - h  qx  sinzot2 

di 


ddq2 

dqx 


PkC, 
An2  a 
dA 


[(-3  sin  A  +  sin  3/1)  Ur  +  (5  cos  A  -  cos  3 A)  Ut  ] 
ddCl 


+q2X - cos  idCl  -  qx  cos  i 

dqx  dqx 


dSq2  PkCBr,  -9  9\f7  /  •  9  •  93W71  dA 

- -  =  — r-2-  (3  cos  A  +  cos3A)Ur  +( sin  A  +  sin  3A)Ut  \  +  qn - ^cosz 

Psi  zlnZ/ if  I— '  '  '  '  —I  Pn 


dSCl 


dq2 

dSq2 

~dQ 


d8q2 

ac7 


An  a 
PkC B 
An2  a 


[(sin  2 A  -  3q2  cos  A  -  3qx  sin  A  -  q2  cos  3 A  +  qx  sin  3  A) 


dq2  dq2 

dO 


dQ 


dU 

+  (-  cos  2 A  +  5qx  cos  A  -  q2  sin  A  -  qx  cos  3 A  -  q2  sin  3 A) 


-qx  cosz 


.dm 

In 


Pk 

An2a 


[(sin  2A  -  3q2  cos  A  -  3qx  sin  A  -  q2  cos  3 A  +  qx  sin  3 A)  Ur 


+ (-  cos  2  A  +  5  qx  cos  A-q2sinA- qx  cos  3  A  -  q2  sin  3A)Ut~^-qx  cos  z 


.dm 

~dC~D 


(C.13) 
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Derivatives  of  m 


_  |_2  sin  A- <7,  sin2/l  +  <7.,cos2/l) 

1  2n2asinr  1  2  ’ 

dm  PkCJJ 


\ 


~  1  1 

- -  .  r1  h  2cos/l  +  —  q,  cos2/t  +  —  q.  sin2/l 

da  n2 a2  sin  A  2  ^  2  j 


dSCl  _  dA 

~W=  lx~de 

dm  PkC  (  -  .  dU  \ r 

di  2n2  asini  \  di  ) 


2cos/l  +  — <7,  cos2/l  +  —qn  sin2/l 
2Hx  2Hi  j 


dm 

dqx 

dm 


PkCBU,  ^  dA 

—  -  -  cos  2/1 +  Q  - — 
An asini 


PkCJJ,  .  ^ 

- - —  B  *  sin2A  +  Q. 

3^  4m  a  sin  z 


n 

dm 

~dn 

dm 


dq, 

dA 

dq2 


PkCn 


2n2  asini 
PkU, 


2cosA  +  —q,cos2A  +  —q^sin2A  — — 
2y'  2Hl  )  dCl 


dCB  2n2  asini 


2 cos  A  +  ^qx  cos2A  +  ^q2  sin2  A 


where 


(C.14) 
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,  P*cn 

x  n2  a 


dSA  2  PkCg  p 


(cos  A  +  2 qx  cos  2  A  +  2 q2  sin  2A}Ur  -  (-  sin  A  +  2 q2  cos  2/1  -  2^  sin  2/1 )  £/  J 
^  2  2  (sinA  +  qx sin  2  A  —  q2  co§2  A^U r  —  (cos/i.  +  q2  sin 2 A  +  qx cos2  A^U t J 


Ml 

da 


-cos  Z 


S3/t  3/1  3<X2 
- =  A, - cosz 

dd  x  dO  d6 


dSA  PkCr: 


di 


n2a 


(sin /l  +  qx  sin 2 A  -  q2  cos2/l)-^L  -(cos/l  +  q2  sin 2 A  +  qx  cos2/l)^7L 


di 


+<X2  sin  z- cosz 


.ddkl 


di 


dSA  PkC  /  -  .  -  x  3/1  .dxi 

- =  — — ^ (U  sin 2A -U tcos2A)  +  A. - cosz - 

n2ay  r  1  ’  1  dqx  dqx 


dq{ 

dSA 

dq2 


PkCB  i  -  -  .  >  dA  .8X1 

-B-\Ur  cos  2  A  +  Ut  sin  2  A  J  +  Ax - cos  z - 


n2  a 


dSA  PkCD 


dd 


n2  a 


dq2 

/  ,dU  i  x3t7 

^sin/l  +  qx  sin 2 A  -  q2  cos2Aj-^-  -  (cos /l  +  q2  sin 2 A  +  qx  cos 2/1 


dq2 
U 

_ r 

dkl 


-cosz- 


dSA  Pk 


.Ml 

3Q 


=  — (sin/l  +  qx  sin2/l  -  q2  cos2A}Ur  -  (cos /l  +  q2  sin2/l  +  qx  cos2/l)L/  J 
dXl 


dCB  n2  a 


dCn 


-cosz 


(C.15) 


Derivatives  of  SCD 

88Cn  d8CR  dSCR  dSCR  88Cn  d8CR 

da  dd  di  dqx  dq2  dCB 

where 


(C.16) 
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dU 

_ r_ 

di 

dU 


=  0 


3Q 

dU 

_ l 

di 

dU. 


=  -sinQcos/l  +cosQsin/l  cos£ 

e  e 

-  =  sin  Q  sin /cos  X  -  cos  Q  sin  /sin  X  cos  £•  +  cos /sin /l  sinf 


dQ 

dUi 

di 

an 


=  -cos  Q  cos  /cos  A  -  sin Q cos  /sin  X  coss 

e  e 

-  =  sin Q  cos  /cos  A  -  cos  Q cos  /sin  A  cos  s-  sin  /sin  A  sin£- 


=  cos  Q  sin /cos  X  +  sin  Q  sin  /sin  X  cos^ 


(C.17) 
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APPENDIX  D.  O  AND  D  MATRIX  FOR  DRAG  PERTURBATIONS 

D.l  0  Matrix 

Derivatives  for  a 


a  =  an  +  a  At 


da 

da0 

da 

de„ 


=  \-anCDCBpp 


‘  3/1  ' 

a—L  +  ~  /, 

V  M  2 


,  dp 

At-a2nCDCBflj^At 


=  0 


=  0 


da 
di0 

-=r~  =  -a2flCDCBpP  If  ^rAt~  a2nCDCB.f^At 

Mo  deo  Mo  Mo 


da  2  ^  ^  „  dfl  M  A.  2  ^  ^  r  dPp  A, 

=  -a  nCDCBp  -L—e-At-a  nC  C  f  —^At 

Mo  M  Mo  M 

da 
3Q 


*20 


=  0 


0 


3a 


dC 


■  -a2nCDppfxAt 


B  0 


(D.l) 


Derivatives  for  i 


—  =  1 

di0 

di  8i  di  di  di  di 

dao  c)0{]  dqw  dq2(]  dQ  0  o'Cfi0 


(D.2) 
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Derivatives  for  q\ 


qx  =  ecos  (O  =  (e0  +  eA^cos£d0  = 

d(h  aM 


^4; 


V 


da( 

dq 


=  <lwAt 


da„ 


L  =  0 


Mo 

lf=0 

a^io  3e0 

M4 


10 


20 


3e0  dq20 


dq 


dn0 

dq. 


L_  _ 


=  0 


3C 


-anCDri2ppf4q  At 


^10  (l"*"-^3^)^10 


0  y 


50 


(D.3) 
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Derivatives  for  qi 


q2  =  esinco  =  (e()  +  eAt)sinft)0  = 

d(h  _  „ 


f  ■  \ 

l  +  —  At 


V 


dan 


^20A^ 


dan 


dq 


30„ 


2  _ 


=  0 


^20  (l  +  /3A?  jt/2o 


"0  J 


dq 


1-0 


di0 

^2 


dq 


=  <hoAt 


10 


de0  dqw 


dq 

dq. 


de0  dq. 


20 


dQ 

dq2 

dC 


^  =  0 


=  -anC  rfp  f  q  At 


BO 


(D.4) 


Derivatives  for  Q 


Q  =  Q0 

dQ  _ 
dQ0 

dQ  _  dQ  _  dQ  _  dQ  _  dQ  _  dQ 

da0  d0{)  di0  dql0  dq20  dCB0 


(D.5) 
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Derivatives  for  0 


where 


30 

dan 


OA  3G  3 q. 

—  A  t- 

\dao 


dG  dq2 
dqt  da0  dq2  daQ 


dO  dG 


30„  30„ 


'dG 

dO 


30  _ 
3/0 

30 

d(ho 

00 

d(ho 

00 

on0 

d6 


OA  At  dG  3gi 
OL  dq  3/ 


V  ""o 
A 


3A 
5<?10 
3  A 


Af- 


A  t- 


^10 

3G 


20 


dq 


10 


=  0 


3C 


=  0 


so 


4  =  [4  +2e/,  +{£'4',,  +  /,)]exp(-/?ae) 

.4 = [4 + H4 + 4) + b2  (34 + 4)>*p(-H 

f 

A  =  ~ankq2pp  =  -ankrf  ppfA 

f4  =  i[(l  +  ap)h  +  (l  -  afi)l2  +  ie(3A  +  ^)]exP(“^ae) 


(D.6) 


(D.7) 


Derivatives  of  /,/3,/4 

%-  =  pe\- jq  +/[  +2(e/0 -Ij/3a-ell)  +  \e{2ell -el0-el2  -2I2/ j3a)]exp(-/3ae) 

=  Mn  (D-8) 

^~  =  aPfn  +[2/i  +f  e(/0  +/2)]exp(-/?«e) 


Approved  for  public  release;  distribution  is  unlimited. 

112 


Ml 

da 


3  _  1 


nCDCBri  ppf^  anCDCB?j  p 


TT  =  1anCDCBep  /4  -  anCDCBqlp  -  anCDCHifj\ 


M L 

p  da 

df4 


■anCDCBTj  fA 


de 


de 


da 

d_Pp 

de 


(D.9) 


^-  =  }[(i-e  +  {e2 -pae)pl0  +  (lp-^al  -j pe)elx  + 

ip  -  2a~l  +  j Pe 2  —  Pe  +  p2ae) I2  - \ (3a_1  +  e/3 ] exp ( ~Pae ) 

V  7  7  (D.10) 

■^  =  y[(-{/fae  +  {/#2a2e  + je-/#a-l)/?a/0  +(2- \e}Palx  Pae  +  \)I2 

+  (y /?ae  -  } p1  a1  e  +  \ e  +  Pa  - 1)  Pal2  -  -jy  (l  ~  Pa) P1 a2 el A ]  exp  (~Pae) 


D.2  D  Matrix 
Derivatives  of  Sa 


dSa 

da 

dSa 

dSa 

di 

dda 

dqx 

dSa 

dq2 

dSa 

In 

dSa 

^C~B 


^  dda,  dSa ,  ddp  ' 


-  +  - 


da  dp  da 

X  ” 

c  dda2  dA 
■  da, — 

1  dA  de 


j 


dSa2 

da2  H - -da, 

2  da  1 


ddat  ddal  ddpp 


dql  dpp  dql 

dda^  ddat  ddp 
-  +  - 


\So,2 


dda2  dda2  dA 


dq2  dp  dq2 


(5*^2  + 


^  dqx  dA  dqxj 

7  dSa2  dda2  dA  ^ 
v  dq2  dA  dq2  y 


da. 


dax 


=  0 


=  -a2CDpp  exp  (~ge)da2 


where 


(Dll) 
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8a  =  8a{8a2 


8a{  =  -a2CDCBpp  exp(-^e) 

8a2=^  +  3^qlsinA-q2  cos  A^j  +  -  tp  (g2  -  q\ )  sin  2/1  - — %2qlq2  cos  2  A 


(D.12) 


88  a 
da 
88a 


~  =  aCDCBPp  eXP(-&)(&-  2) 


dqx 

88a 


8q 

88a 


1  =  a1C[)CBpp8  exp  (-^e)  cos  co 
1  =  a  2CDCB/?^  exp  (- sin  to 


dp 


-  =  -a2cDCB  exp(-^e) 


<3c>a2  1 

da  ~  H 
88 a 


™PP 

Po  ( 

=  — -exp  - 

1 

1 

1 

da 

H  \ 

h  ;v  7 

8SPp 

f 

r-r '  ^ 

8q{ 

=  /V?exP^- 

— -  costy 

H  J 

88p 

r 

rn~r(^ 

dq2 

=  /V?exP^_ 

— -  sm<y 

H  J 

A-q2 

COS  /l)  +  ^ 

’  4  H 

(c/2  -q2jsin2A 

_ 

2h' 


8A 

8A 


- +  3^j  cos/L  +  q2 sin/l) -\ - [q2  - g2 )  cos  2/1 -\ - qtq2 sin 2 A 

1 


86  1  +  2^  cos  A  +  2q2  sin  A 

H=X- 

P 


(D.13) 


(D.14) 


Derivatives  of  8i 


88i  _  88i  88i  _  88i  88i  _  88i  88i 
da  86  8i  8qx  dq2  8  Q  8CB 


(D.15) 


Approved  for  public  release;  distribution  is  unlimited. 

114 


Derivatives  of  8qx 


,  dd3n  dSpp\3a 

da  J  1 


da 

dSq^ 

dd 

dSqx 


=  8qx 


=  0 


da  dpp 
dSqx2  dA 

dl  dP 


+ 


^12 

da 


8q. 


di 

dSqx  ( dSqx 


dqx  {  dqx 

d_Sq_  i  f  dSq , 
dq2 

dSq^ 

d  Q 


+ 


dq 

dq2  dpD  dq2  ) 


-  +  - 


+ 


f  ^12  ,  ^12  5/1  1 


dq2 


+ 


ds l 


dqj5qi 


=  0 


=  -\aCDppexp(-Ze)SqX2 


(D.16) 


where 


Sqx  =  8qn8qx2 

8qxx  =  -  ^  aCDCBpp  exp(-^e) 

<fy2  =  ^2  + 1 £  V  +  ^ ^2)  sinA  ~  \ cos  A  + 

sin  2/1  -  q2  cos  2/l)  +  ^  (g2  -  q2  j  sin  3 A  -  ^  <^2qtq2  cos  3/1 


(D.17) 


dggn 

5a 

5^! 

dq2 

d8qn 

dpp 


~\CdCbPpq* p(-fr)(l-fr) 

\aCDCBpp%Qx  p(-£e)  COS  (2? 

2aC«C«^/exP(“^)  sin  <2? 
-^CZ)C5eXP(-^) 


(D.18) 
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dSq, 


12 


da 


“U*I+2*2 

I_ 

2  H 


—  sin  2  -  —q.q^  cos  2  + 

H  HHxHl 

X—  [qx  sin 2/1  -  q2  cos 2 A)  +  ~z{q\  ~  q] )  sin  32  -  ^  ^  Wi cos  3A 


dSq, 


12 


=  2  +  j€24?+-7?tf  cosA  +  -^qiq2  sin 2  + 


-'2  2 


dA 


dSq 


12 


dqx 


dSq , 


12 


dq2 


(<^+ 3  cos  22  +  q2  sin  22)  +  (<?2  _  ^  )cos32  +  sin  32 

=  —(?>qx  sin  2  -  g2  cos  2)  +  ^  sin22  +  ~~(qt  sin  32  -  q2  cos  32) 

&+3  &2 

=  — (g2  sin  2  -  ^  cos  2) - cos  22 - sin  32  +  qx  cos  32) 


CD-19) 
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dSn  dSD.  dSD.  dSD.  dSD.  dSD.  dSD. 

da  d6  di  dq]  dq2  d  Q  dCB 


(D.20) 


(D.21) 
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Derivatives  of  86 


886 

da 

886 

86 


886_ 

8i 

886 

8qx 

886 

dq2 

886 

80 

886 

dCB 


^1,0  0  ,  O  •  o  •  a  dSch  i  dSch 

=  1 1  +  2g  cos  8  +  2q  2  sin  8) +2sin/t -  2cos/l - 

v  '  da  da  8a 

/'i  ,  t  o  ,  o  •  i\ss^  o  •  ,  dS(h  o  ,  88q2 

■1  +  2(7,  cos  A  +  2q^  sin  A] +2sm/t -  2cos/l - 

1  ^2  ;  a#  ^  dn 


_/  .  ,\3/l  _  .5/ 1  _  _  .  .  5/ 1  _ 

+2  -o,  sin/l  +  o^  cos /l  — o/t+2cos/l — bq ,  +2sin/l  —  do,, 
V  *i  ^2  3#  yi  50  y2 

li  ,  o  o  ,  o  •  o  •  o  ,  88 q 2 

=  1  +  2  o,  cos  A  +  2o  sin  A) - +2sin/t - -  -2cos8 - - 

v  ’da  8i  8i 

(.  .  ...  3d/l  .  .  , 

=  1  +  2g  cos  X  +  2q  sm  X  - +2smX - 

v  8q  8q 


1  -  2cos/l  — — +  28X  cos  2 


8qx 


-  (l  +  2  q,  cos /l +  2  q.  smA\^^+2smX<—^  -  2cos2  +  28X  sin  A 

1  y2  )  8  fs 


dq2 

80. 


=  (l  +  2o,  cos /l +  2 q,  sin  8]^^L+28m8  — — —  -  2cos2 
1  y2  ;  aQ  aQ 

li  ,  o  ,  ,  o  •  |0  .  ,  88q  88 q 

=  1  +  2 o,  cos  A  +  2o,,  sm  A  +2sm/l L  -  2cos A 

1  1  2  'SC„  dC„  5C„ 


(D.22) 


where 


&&  _  f  88Xx  88X1  8Sp)  d8X2 

da  ^  da  8pp  da  J  da 

88X  _  S1  dSX2  88 
86  ~  A'  8X  86 


88X  f  88X.  88X.  88p  )  ( Q8X2  8SX2  5/0  „ 

8qx  y  8qx  8pp  8qx  J  y  8qx  8X  8qx) 

888  (d8X.  88Xx88p  )  f  88X2  8882  88 )  „ 

dq2  l  dq2  +  dPP  d(h  J  2+l  dch  +  dA  d(h )  1 


888 

80 

888 

scB 


=  0 

=  ^aCDPpeM~^e)(U2 


(D.23) 


and 
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(D.24) 


8A  =  8AX8A2 

S^=\aCrfBppQM~&) 

8A2  =  qi  cos  A  +  q2  sin  A 


3SA[ 

da 

dSAl 

dqx 

dSAl 

dq2 

dSAl 

dpp 


\CdCbPpq*p{-&)(\-&) 

-\aCDCBPP^QXP{-^)cos(o 

~\aCDCBpp^ p(-£e)sin© 
\uCdCbgl p(-£e) 


dSA2 

da 

dSA2 

8A 

d8A2 

dqx 

dSA^ 

dq2 


0 

-ql  sinA  +  q2  cos  A 
cos  A 

sin  A 


(D.25) 


(D.26) 
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LIST  OF  ACRONYMS,  ABBREVIATIONS,  AND  SYMBOLS 


Acronym/  Abbreviation  Description 


AMR 

GA-STM 

GMAT 

GNC 

GUI 

LEO 

LVLH 

MMS 

NASA 

PCO 

RMS 

SRP 

STM 


Area-to  mass  ratio 

Gim  Alfriend-State  Transition  Matrix 

Goddard  General  Mission  Analysis  Tool 

Guidance,  navigation,  and  control 

Graphical  User  Interface 

Low  earth  orbit 

Local  vertical  local  horizontal 

Magnetospheric  Multiscale 

National  Aeronautics  and  Space  Administration 

Projected  Circular  Orbit 

Root  mean  square 

Solar  radiation  pressure 

State  Transition  Matrix 
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